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

Authors: Salvatore Brischetto
Journal: Shock and Vibration
 (2014)

---

## Abstract

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 vibration modes for a number of thickness ratios, imposed wave numbers, geometries, embedded materials, 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).


---

## Body

## 1. 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 three-dimensional 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 angle-ply 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 two-dimensional 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–35]. Exact closed form solutions (based on three-dimensional elasticity theory) were carried out in [33] for both in-plane 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 three-dimensional 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–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.


## 2. 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 thethickness of the shell at that point [42, 43]. Geometry and the reference system are indicated in Figure 1 where curvilinear orthogonal coordinates (α, β, z) are shown. The three-dimensional space (the space in which the three independent variables α, β, and z vary) is an Euclidean space. Rα and Rβ are the principal radii of curvature along the coordinates α and β, respectively. z is the thickness coordinate that varies from -h/2 to +h/2 (where h is the thickness of the structure as shown in Figure 1).Figure 1
Notation and reference system for shells.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–46]):(1)ϵαα=1Hα∂u∂α+wHαRα,ϵββ=1Hβ∂v∂β+wHβRβ,ϵzz=∂w∂z,γαβ=1Hα∂v∂α+1Hβ∂u∂β,γαz=1Hα∂w∂α+∂u∂z-uHαRα,γβz=1Hβ∂w∂β+∂v∂z-vHβRβ,
where u, v, and w are the displacement components in α, β, and z 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
(2)Hα=(1+zRα),Hβ=(1+zRβ),Hz=1.
General geometrical relations for spherical shells degenerate into geometrical relations for cylindrical shells when Rα or Rβ is infinite (with Hα or Hβ equals one), and they degenerate into geometrical relations for plates when both Rα and Rβ are infinite (with Hα = Hβ = 1) [46].Geometrical relations (see (1)) are substituted in three-dimensional linear elastic constitutive equations in orthogonal curvilinear coordinates (α, β, z) for orthotropic material in the structural reference system (see [47, 48]). Partial derivatives ∂/∂α, ∂/∂β, and ∂/∂z are indicated with subscripts ,α, ,β, and ,z. A closed form solution of differential equations of equilibrium for shells is obtained when coefficients C16, C26, C36, and C45 are set to zero (this means orthotropic angle θ equals 0∘ or 90∘):
(3)σαα=C11Hαu,α+C11HαRαw+C12Hβv,β+C12HβRβw+C13w,z,σββ=C12Hαu,α+C12HαRαw+C22Hβv,β+C22HβRβw+C23w,z,σzz=C13Hαu,α+C13HαRαw+C23Hβv,β+C23HβRβw+C33w,z,σβz=C44Hβw,β+C44v,z-C44HβRβv,σαz=C55Hαw,α+C55u,z-C55HαRαu,σαβ=C66Hαv,α+C66Hβu,β.


## 3. Equilibrium Equations

The three differential equations of equilibrium written for the case of free vibration analysis of spherical shells with constant radii of curvatureRα and Rβ are here given (the most general form for variable radii of curvature can be found in [45, 49]):
(4)Hβ∂σαα∂α+Hα∂σαβ∂β+HαHβ∂σαz∂z+(2HβRα+HαRβ)σαz=ρHαHβu¨,Hβ∂σαβ∂α+Hα∂σββ∂β+HαHβ∂σβz∂z+(2HαRβ+HβRα)σβz=ρHαHβv¨,Hβ∂σαz∂α+Hα∂σβz∂β+HαHβ∂σzz∂z-HβRασαα-HαRβσββ+(HβRα+HαRβ)σzz=ρHαHβw¨,
where ρ is the mass density and u¨, v¨, and w¨ 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 C16=C26=C36=C45=0). The three displacement components have the following harmonic form:
(5)u(α,β,z,t)=U(z)eiωtcos⁡(α¯α)sin(β¯β),v(α,β,z,t)=V(z)eiωtsin(α¯α)cos⁡(β¯β),w(α,β,z,t)=W(z)eiωtsin(α¯α)sin(β¯β),
where U, V, and W are the displacement amplitudes in α, β, and z directions, respectively. i is the coefficient of the imaginary unit, ω=2πf is the circular frequency where f is the frequency value, and t is the time. In coefficients α-=mπ/a and β-=nπ/b, m and n are the half-wave numbers and a and b 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 Rα or Rβ is infinite (with Hα or Hβ equals one) and into equilibrium equations for plates when Rα and Rβ are infinite (with Hα and Hβ equal one). In this way a unique and general closed formulation is possible for any geometry:
(6)(-C55HβHαRα2-C55RαRβ-α¯2C11HβHα-β¯2C66HαHβ+ρHαHβω2)U+(-α¯β¯C12-α¯β¯C66)V+(α¯C11HβHαRα+α¯C12Rβ+α¯C55HβHαRα+α¯C55Rβ)W+(C55HβRα+C55HαRβ)U,z+(α¯C13Hβ+α¯C55Hβ)W,z+(C55HαHβ)U,zz=0,(-α¯β¯C66-α¯β¯C12)U+(-C44HαHβRβ2-C44RαRβi-α¯2C66HβHα-β¯2C22HαHβ+ρHαHβω2C44HαHβRβ2)V+(β¯C44HαHβRβ+β¯C44Rα+β¯C22HαHβRβ+β¯C12Rα)W+(C44HαRβ+C44HβRα)V,z+(β¯C44Hα+β¯C23Hα)W,z+(C44HαHβ)V,zz=0,(α¯C55HβHαRα-α¯C13Rβ+α¯C11HβHαRα+α¯C12Rβ)U+(β¯C44HαHβRβ-β¯C23Rα+β¯C22HαHβRβ+β¯C12Rα)V+(C13RαRβ+C23RαRβ-C11HβHαRα2-2C12RαRβ-C22HαHβRβ2-α¯2C55HβHα-β¯2C44HαHβ+ρHαHβω2C11HβHαRα2)W+(-α¯C55Hβ-α¯C13Hβ)U,z+(-β¯C44Hα-β¯C23Hα)V,z+(C33HβRα+C33HαRβ)W,z+(C33HαHβ)W,zz=0.
The system of equation (6) can be written in a compact form by introducing coefficients Ar for each block ():
(7)A1U+A2V+A3W+A4U,z+A5W,z+A6U,zz=0,A7U+A8V+A9W+A10V,z+A11W,z+A12V,zz=0,A13U+A14V+A15W+A16U,z+A17V,z+A18W,z+A19W,zz=0.
Equation (7) is a system of three second order differential equations. They are written for spherical shell panels with constant radii of curvature but they automatically degenerate into equations for cylindrical shells and plates (see also [7, 25]).


### 3.1. 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:
(8)[A6000000A12000000A19000000A6000000A12000000A19][UVWU′V′W′]′=[000A6000000A12000000A19-A1-A2-A3-A40-A5-A7-A8-A90-A10-A11-A13-A14-A15-A16-A17-A18][UVWU′V′W′].
Equation (8) can be written in a compact form as
(9)D∂U∂z~=AU,
where ∂U/∂z~=U′ and U=[UVWU′V′W′]. Equation (9) can be written as
(10)DU′=AU,(11)U′=D-1AU,(12)U′=A*U,
with A*=D-1A.In the case of plate geometry coefficientsA3, A4, A9, A10, A13, A14, and A18 are zero because the radii of curvature Rα and Rβ are infinite. The other coefficients A1, A2, A5, A6, A7, A8, A11, A12, A15, A16, A17, and A19 are constant because parametric coefficients Hα=Hβ=1 and they do not depend on the thickness coordinate z~; therefore matrices D, A, and A* are constant in the plate case. The solution of (12) can be written as [51, 52]
(13)U(z~)=exp⁡(A*z~)U(0)withz~∈[0,h],
where z~ is the thickness coordinate referred to the bottom of the structure that has values from 0 at the bottom to h at the top (see Figure 2). The exponential matrix is calculated with z~=h as
(14)A**=exp⁡(A*h)=I+A*h+A*22!h2A**=exp(A*h)=+A*33!h3+⋯+A*NN!hN,
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 one-layered 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:
(15)T=I,
Equation (15) allows the calculation of the matrix H as
(16)H=A**I=A**;
this last matrix permits linking vector U calculated at the top of the plate (z~=h) with vector U calculated at the bottom of the plate (z~=0):
(17)U(h)=HU(0).Thickness coordinates for plates and shells.
(a)
(b)In the case of shell geometry matricesD, A, and A* are not constant because of parametric coefficients Hα=(1+z/Rα)=(1+(z~-h/2)/Rα) and Hβ=(1+z/Rβ)=(1+(z~-h/2)/Rβ) that depend on z or z~ coordinate (see Figure 2). A first method could be the use of hypothesis z/Rα=z/Rβ=0 (it is valid only for very thin shells) that means Hα=Hβ=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 Hα=Hβ=1 will be given at the end of the section for results. The second method (used in this paper) is the introduction of several j fictitious layers where Hα and Hβ can exactly be calculated with the z~ thickness coordinate of the middle of the fictitious layers; in this way (12), (13), and (14) are valid for each j layer because Aj* is here constant. The exponential matrix Aj** is calculated in each fictitious j layer with the thickness hj of the layer. The results proposed in this paper for the benchmarks are given with a number of fictitious layers j=M=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 M layers for shell geometry, it is necessary to calculate M-1 transfer matrices Tj-1,j by using for each interface the following conditions:
(18)ujb=uj-1t,vjb=vj-1t,wjb=wj-1t,σzzjb=σzzj-1t,σαzjb=σαzj-1t,σβzjb=σβzj-1t;
that means each displacement and transverse stress component at the top (t) of the j-1 layer is equal to displacement and transverse stress component at the bottom (b) of the j layer. The calculated Tj-1,j matrices allow linking U at the bottom of the j layer with U at the top of the j-1 layer:
(19)Uj(0)=Tj-1,jUj-1(hj-1).U at the top of the j layer is linked with U at the bottom of the same layer j by means of the exponential matrix Aj**:
(20)Uj(hj)=Aj**Uj(0).
Equation (19) can recursively be introduced in (20) for the M-1 interfaces to obtain
(21)UM(hM)TM-1,MAM-1**TM-2,M-1⋯A2**T1,2A1**U1(0)=AM**TM-1,MAM-1**TM-2,M-1⋯A2**T1,2A1**U1(0),
by defining the matrix Hm for the multi-fictitious-layer shells, it is possible to write (21) in analogy with (17):
(22)UM(hM)=HmU1(0).
It links U calculated at the top of the last M fictitious layer with U calculated at the bottom of the first layer. Matrices Aj** are constant because they are evaluated with Rα, Rβ, α¯, and β¯ calculated in the midsurface Ω0 of the shell and with Hα and Hβ calculated in the middle of each j layer. Matrices Tj-1,j are constant because they are evaluated with Rα, Rβ, α¯, and β¯ calculated in the midsurface Ω0 of the shell and with Hα and Hβ calculated at each fictitious interface. Equation (17) is valid for one-layered plates and (22) is valid for one-layered shells divided in j 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:(23)σzz=σαz=σβz=0,forz=-h2,+h2orz~=0,h,(24)w=v=0,σαα=0,forα=0,a,(25)w=u=0,σββ=0,forβ=0,b.
Equation (23) imposed at the top of the structure has Rα, Rβ, α¯, and β¯ calculated in the midsurface Ω0 of the shell, and Hαt and Hβt are calculated at the top of the shell for z~=h:
(26)σzzt=-α¯C13HαtUt+C13HαtRαWt-β¯C23HβtVt+C23HβtRβWt+C33W,zt=0,σβzt=β¯C44HβtWt+C44V,zt-C44HβtRβVt=0,σαzt=α¯C55HαtWt+C55U,zt-C55HαtRαUt=0.
Equation (23) imposed at the bottom of the structure has Rα, Rβ, α¯, and β¯ calculated in the midsurface Ω0 of the shell, and Hαb and Hβb are calculated at the bottom of the shell for z~=0:
(27)σzzb=-α¯C13HαbUb+C13HαbRαWb-β¯C23HβbVb+C23HβbRβWb+C33W,zb=0,σβzb=β¯C44HβbWb+C44V,zb-C44HβbRβVb=0,σαzb=α¯C55HαbWb+C55U,zb-C55HαbRαUb=0.
Equations (26) in matrix form are (UM(hM) means U calculated at the top of the shell)
(28)[-α¯C13Hαt-β¯C23Hβt(C13HαtRα+C23HβtRβ)00C330-C44HβtRββ¯C44Hβt0C440-C55HαtRα0α¯C55HαtC5500]×[UM(hM)VM(hM)WM(hM)UM′(hM)VM′(hM)WM′(hM)]=[000].
Equation (27) in matrix form is (U1(0) means U calculated at the bottom of the shell)
(29)[-α¯C13Hαb-β¯C23Hβb(C13HαbRα+C23HβbRβ)00C330-C44HβbRββ¯C44Hβb0C440-C55HαbRα0α¯C55HαbC5500]×[U1(0)V1(0)W1(0)U1′(0)V1′(0)W1′(0)]=[000].
Equations (28)-(29) in compact form to express the free stress state at the top and bottom of the shell are
(30)BM(hM)UM(hM)=0,(31)B1(0)U1(0)=0;
in the case of plate, there are no j fictitious layers:
(32)B(h)U(h)=0,B(0)U(0)=0.
Equation (22) can be substituted in (30):
(33)BM(hM)HmU1(0)=0;B1(0)U1(0)=0,
in this way (33) are grouped in the following system:
(34)[BM(hM)HmB1(0)][U1(0)]=[0],
and introducing the (6×6) E matrix, (34) becomes:
(35)[E][U1(0)]=[0].
Equation (35) is also valid for plate case where the fictitious layers are not introduced and U(0)=U1(0), H=Hm, and B=BM(hM)=B1(0). The solution is implemented in a Matlab code where only the spherical shell method is considered; it automatically degenerates into cylindrical open/closed shell and plate methods.The free vibration problem means finding the nontrivial solution ofU1(0) in (35); this means imposing the determinant of matrix E equal to zero:
(36)det⁡[E]=0.
Equation (36) needs to find the roots of a higher order polynomial in λ=ω2. For each pair of half-wave numbers (m,n) a certain number of circular frequencies are obtained depending on the order N chosen for the exponential matrices Aj** and the M number of fictitious layers.


### 3.2. Vibration Modes

A certain number of circular frequenciesωl are found when half-wave numbers m and n are imposed in the structures. For each frequency ωl, it is possible to find the vibration mode through the thickness in terms of the three displacement components. If the frequency ωl 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 ωl, the vector U calculated at the bottom:
(37)U1ωl(0)=[U1(0)V1(0)W1(0)U1′(0)V1′(0)W1′(0)]ωlT.T means the transpose of the vector and the subscript ωl means that the null space is calculated for the circular frequency ωl.It is possible to findUjωl(z~j) (with the three displacement components Ujωl(z~j), Vjωl(z~j), and Wjωl(z~j) through the thickness) for each j layer by using (19)–(22); the thickness coordinate z~ 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 j fictitious layers and the M number equals 1.


## 3.1. 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:
(8)[A6000000A12000000A19000000A6000000A12000000A19][UVWU′V′W′]′=[000A6000000A12000000A19-A1-A2-A3-A40-A5-A7-A8-A90-A10-A11-A13-A14-A15-A16-A17-A18][UVWU′V′W′].
Equation (8) can be written in a compact form as
(9)D∂U∂z~=AU,
where ∂U/∂z~=U′ and U=[UVWU′V′W′]. Equation (9) can be written as
(10)DU′=AU,(11)U′=D-1AU,(12)U′=A*U,
with A*=D-1A.In the case of plate geometry coefficientsA3, A4, A9, A10, A13, A14, and A18 are zero because the radii of curvature Rα and Rβ are infinite. The other coefficients A1, A2, A5, A6, A7, A8, A11, A12, A15, A16, A17, and A19 are constant because parametric coefficients Hα=Hβ=1 and they do not depend on the thickness coordinate z~; therefore matrices D, A, and A* are constant in the plate case. The solution of (12) can be written as [51, 52]
(13)U(z~)=exp⁡(A*z~)U(0)withz~∈[0,h],
where z~ is the thickness coordinate referred to the bottom of the structure that has values from 0 at the bottom to h at the top (see Figure 2). The exponential matrix is calculated with z~=h as
(14)A**=exp⁡(A*h)=I+A*h+A*22!h2A**=exp(A*h)=+A*33!h3+⋯+A*NN!hN,
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 one-layered 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:
(15)T=I,
Equation (15) allows the calculation of the matrix H as
(16)H=A**I=A**;
this last matrix permits linking vector U calculated at the top of the plate (z~=h) with vector U calculated at the bottom of the plate (z~=0):
(17)U(h)=HU(0).Thickness coordinates for plates and shells.
(a)
(b)In the case of shell geometry matricesD, A, and A* are not constant because of parametric coefficients Hα=(1+z/Rα)=(1+(z~-h/2)/Rα) and Hβ=(1+z/Rβ)=(1+(z~-h/2)/Rβ) that depend on z or z~ coordinate (see Figure 2). A first method could be the use of hypothesis z/Rα=z/Rβ=0 (it is valid only for very thin shells) that means Hα=Hβ=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 Hα=Hβ=1 will be given at the end of the section for results. The second method (used in this paper) is the introduction of several j fictitious layers where Hα and Hβ can exactly be calculated with the z~ thickness coordinate of the middle of the fictitious layers; in this way (12), (13), and (14) are valid for each j layer because Aj* is here constant. The exponential matrix Aj** is calculated in each fictitious j layer with the thickness hj of the layer. The results proposed in this paper for the benchmarks are given with a number of fictitious layers j=M=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 M layers for shell geometry, it is necessary to calculate M-1 transfer matrices Tj-1,j by using for each interface the following conditions:
(18)ujb=uj-1t,vjb=vj-1t,wjb=wj-1t,σzzjb=σzzj-1t,σαzjb=σαzj-1t,σβzjb=σβzj-1t;
that means each displacement and transverse stress component at the top (t) of the j-1 layer is equal to displacement and transverse stress component at the bottom (b) of the j layer. The calculated Tj-1,j matrices allow linking U at the bottom of the j layer with U at the top of the j-1 layer:
(19)Uj(0)=Tj-1,jUj-1(hj-1).U at the top of the j layer is linked with U at the bottom of the same layer j by means of the exponential matrix Aj**:
(20)Uj(hj)=Aj**Uj(0).
Equation (19) can recursively be introduced in (20) for the M-1 interfaces to obtain
(21)UM(hM)TM-1,MAM-1**TM-2,M-1⋯A2**T1,2A1**U1(0)=AM**TM-1,MAM-1**TM-2,M-1⋯A2**T1,2A1**U1(0),
by defining the matrix Hm for the multi-fictitious-layer shells, it is possible to write (21) in analogy with (17):
(22)UM(hM)=HmU1(0).
It links U calculated at the top of the last M fictitious layer with U calculated at the bottom of the first layer. Matrices Aj** are constant because they are evaluated with Rα, Rβ, α¯, and β¯ calculated in the midsurface Ω0 of the shell and with Hα and Hβ calculated in the middle of each j layer. Matrices Tj-1,j are constant because they are evaluated with Rα, Rβ, α¯, and β¯ calculated in the midsurface Ω0 of the shell and with Hα and Hβ calculated at each fictitious interface. Equation (17) is valid for one-layered plates and (22) is valid for one-layered shells divided in j 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:(23)σzz=σαz=σβz=0,forz=-h2,+h2orz~=0,h,(24)w=v=0,σαα=0,forα=0,a,(25)w=u=0,σββ=0,forβ=0,b.
Equation (23) imposed at the top of the structure has Rα, Rβ, α¯, and β¯ calculated in the midsurface Ω0 of the shell, and Hαt and Hβt are calculated at the top of the shell for z~=h:
(26)σzzt=-α¯C13HαtUt+C13HαtRαWt-β¯C23HβtVt+C23HβtRβWt+C33W,zt=0,σβzt=β¯C44HβtWt+C44V,zt-C44HβtRβVt=0,σαzt=α¯C55HαtWt+C55U,zt-C55HαtRαUt=0.
Equation (23) imposed at the bottom of the structure has Rα, Rβ, α¯, and β¯ calculated in the midsurface Ω0 of the shell, and Hαb and Hβb are calculated at the bottom of the shell for z~=0:
(27)σzzb=-α¯C13HαbUb+C13HαbRαWb-β¯C23HβbVb+C23HβbRβWb+C33W,zb=0,σβzb=β¯C44HβbWb+C44V,zb-C44HβbRβVb=0,σαzb=α¯C55HαbWb+C55U,zb-C55HαbRαUb=0.
Equations (26) in matrix form are (UM(hM) means U calculated at the top of the shell)
(28)[-α¯C13Hαt-β¯C23Hβt(C13HαtRα+C23HβtRβ)00C330-C44HβtRββ¯C44Hβt0C440-C55HαtRα0α¯C55HαtC5500]×[UM(hM)VM(hM)WM(hM)UM′(hM)VM′(hM)WM′(hM)]=[000].
Equation (27) in matrix form is (U1(0) means U calculated at the bottom of the shell)
(29)[-α¯C13Hαb-β¯C23Hβb(C13HαbRα+C23HβbRβ)00C330-C44HβbRββ¯C44Hβb0C440-C55HαbRα0α¯C55HαbC5500]×[U1(0)V1(0)W1(0)U1′(0)V1′(0)W1′(0)]=[000].
Equations (28)-(29) in compact form to express the free stress state at the top and bottom of the shell are
(30)BM(hM)UM(hM)=0,(31)B1(0)U1(0)=0;
in the case of plate, there are no j fictitious layers:
(32)B(h)U(h)=0,B(0)U(0)=0.
Equation (22) can be substituted in (30):
(33)BM(hM)HmU1(0)=0;B1(0)U1(0)=0,
in this way (33) are grouped in the following system:
(34)[BM(hM)HmB1(0)][U1(0)]=[0],
and introducing the (6×6) E matrix, (34) becomes:
(35)[E][U1(0)]=[0].
Equation (35) is also valid for plate case where the fictitious layers are not introduced and U(0)=U1(0), H=Hm, and B=BM(hM)=B1(0). The solution is implemented in a Matlab code where only the spherical shell method is considered; it automatically degenerates into cylindrical open/closed shell and plate methods.The free vibration problem means finding the nontrivial solution ofU1(0) in (35); this means imposing the determinant of matrix E equal to zero:
(36)det⁡[E]=0.
Equation (36) needs to find the roots of a higher order polynomial in λ=ω2. For each pair of half-wave numbers (m,n) a certain number of circular frequencies are obtained depending on the order N chosen for the exponential matrices Aj** and the M number of fictitious layers.


## 3.2. Vibration Modes

A certain number of circular frequenciesωl are found when half-wave numbers m and n are imposed in the structures. For each frequency ωl, it is possible to find the vibration mode through the thickness in terms of the three displacement components. If the frequency ωl 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 ωl, the vector U calculated at the bottom:
(37)U1ωl(0)=[U1(0)V1(0)W1(0)U1′(0)V1′(0)W1′(0)]ωlT.T means the transpose of the vector and the subscript ωl means that the null space is calculated for the circular frequency ωl.It is possible to findUjωl(z~j) (with the three displacement components Ujωl(z~j), Vjωl(z~j), and Wjωl(z~j) through the thickness) for each j layer by using (19)–(22); the thickness coordinate z~ 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 j fictitious layers and the M number equals 1.


## 4. 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 one-layered 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 Hα and Hβ parametric coefficients have also been investigated.Geometries considered for the assessments and benchmarks: (a) square plate, (b) rectangular plate, (c) cylindrical shell panel, (d) cylinder, and (e) spherical shell panel.
(a)
(b)
(c)
(d)
(e)


### 4.1. 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 E=70GPa, Poisson ratio ν=0.3, and mass density ρ=2702kg/m3. The thickness value is h=1. Two different thickness ratios are considered, that is, a/h=10 and a/h=10. The first six circular frequencies in no-dimensional form ω-=ω(a2/h)ρ/E are given in Table 1 for imposed half-wave numbers m=n=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 M number of fictitious layers and the N 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 Hα and Hβ. For M=1, an order of expansion N=18 for the exponential matrix always gives the correct solution for each vibration mode and thickness ratio a/h. The calculation of the exponential matrix is not time consuming. An alternative method could be the use of M fictitious layers and a lower order of expansion N=3 for the exponential matrix. In this last case, M=100 fictitious layers always give the correct solution for each vibration mode and thickness ratio. The use of M=100 layers is not time consuming because the E matrix has always 6×6 dimension for any value of M layers, even if the proposed method uses a layerwise approach. Solutions with M=1 and N=18 or with M=100 and N=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.Table 1
Assessment 1. Simply supported aluminum square plate with Poisson ratioν=0.3. First six exact natural circular frequencies in no-dimensional form ω-=ω(a2/h)ρ/E for imposed half-wave numbers m=n=1. Comparison between present three-dimensional analysis and three-dimensional analyses by Vel and Batra [34] and Srinivas et al. [22]. M is the number of fictitious layers and N is the order of expansion for the exponential matrix.
Mode
I
II
III
IV
V
VI
a
/
h
=
10
Present 3D (M=1, N=3)
5.3370
27.554
46.493
154.39
168.09
268.36
Present 3D (M=1, N=7)
5.7769
27.554
46.502
192.91
197.98
329.35
Present 3D (M=1, N=10)
5.7769
27.554
46.502
196.21
201.87
308.34
Present 3D (M=1, N=15)
5.7769
27.554
46.502
196.77
201.34
358.23
Present 3D (M=1, N=18)
5.7769
27.554
46.502
196.77
201.34
357.42
Present 3D (M=10, N=3)
5.7769
27.554
46.502
196.71
201.27
355.40
Present 3D (M=50, N=3)
5.7769
27.554
46.502
196.77
201.34
351.90
Present 3D (M=100, N=3)
5.7769
27.554
46.502
196.77
201.34
357.42
3D [34]
5.7769
27.554
46.503
196.77
201.34
357.42
3D [22]
5.7769
27.554
46.502
196.77
201.34
357.42
a
/
h
=
10
Present 3D (M=1, N=3)
/
8.7132
14.182
17.512
/
/
Present 3D (M=1, N=7)
4.6598
8.7132
14.463
20.987
23.999
29.925
Present 3D (M=1, N=10)
4.6582
8.7132
14.463
21.383
25.009
31.399
Present 3D (M=1, N=15)
4.6581
8.7132
14.463
21.343
24.830
34.049
Present 3D (M=1, N=18)
4.6581
8.7132
14.463
21.343
24.830
33.982
Present 3D (M=10, N=3)
4.6580
8.7132
14.463
21.337
24.821
33.939
Present 3D (M=50, N=3)
4.6581
8.7132
14.463
21.343
24.830
33.923
Present 3D (M=100, N=3)
4.6581
8.7132
14.463
21.343
24.830
33.982
3D [34]
4.6582
8.7132
14.463
21.343
24.830
33.982
3D [22]
4.6582
8.7134
14.462
21.343
24.830
33.982The 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 Rα=1 and infinite radius of curvature Rβ in β direction. The circular dimension is a=2πRα; the thickness h is 0.12 and 0.18; that means thickness ratios are h/Rα=0.12 and h/Rα=0.18, respectively. Several ratios cylinder length/radius of curvature (b/Rα) are investigated in the case of b=2,1,0.5,0.25. The fundamental circular frequency in no-dimensional form ω-=ω(h/π)ρ/G is given in Tables 2 and 3 (thickness ratio h/Rα equals 0.12 ad 0.18, resp.) for axial half-wave number n=1 and several circular half-wave numbers m=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 Hα and Hβ. The convergence study is shown for a value N=3 for the expansion of the exponential matrix. M=100 fictitious layers always give the correct solution for each thickness ratio h/Rα (see Tables 2 and 3), length b of the cylinder, and half-wave numbers m in circumferential direction. The present method is not time consuming; M=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.Table 2
Assessment 2 (part I). Simply supported aluminum cylinder with Poisson ratioν=0.3 and thickness ratio h/Rα=0.12. Lowest exact natural circular frequencies in no-dimensional form ω-=ω(h/π)ρ/G for imposed axial half-wave number n=1 and several half-wave numbers m in circumferential direction. Comparison between present three-dimensional analysis and three-dimensional analysis by Armenakas et al. [8] as adapted in Bhimaraddi [9]. M is the number of fictitious layers and N is the order of expansion for the exponential matrix.
h
/
R
α
=
0.12
m
=
2
m
=
4
m
=
6
m
=
8
b
/
R
α
=
2
Present 3D (M=1, N=3)
0.03724
0.02323
0.02175
0.02846
Present 3D (M=10, N=3)
0.03730
0.02360
0.02464
0.03687
Present 3D (M=50, N=3)
0.03730
0.02359
0.02463
0.03686
Present 3D (M=100, N=3)
0.03730
0.02359
0.02462
0.03686
3D [8, 9]
0.03730
0.02359
0.02462
0.03686
b
/
R
α
=
1
Present 3D (M=1, N=3)
0.05769
0.04772
0.04229
0.04272
Present 3D (M=10, N=3)
0.05853
0.04979
0.04790
0.05546
Present 3D (M=50, N=3)
0.05853
0.04978
0.04789
0.05545
Present 3D (M=100, N=3)
0.05853
0.04978
0.04789
0.05545
3D [8, 9]
0.05853
0.04978
0.04789
0.05545
b
/
R
α
=
0.5
Present 3D (M=1, N=3)
0.08017
0.07732
0.07308
0.06667
Present 3D (M=10, N=3)
0.10057
0.10234
0.10688
0.11528
Present 3D (M=50, N=3)
0.10057
0.10234
0.10688
0.11528
Present 3D (M=100, N=3)
0.10057
0.10234
0.10688
0.11528
3D [8, 9]
0.10057
0.10234
0.10688
0.11528
b
/
R
α
=
0.25
Present 3D (M=1, N=3)
0.48156
0.48619
0.49378
0.50419
Present 3D (M=10, N=3)
0.27489
0.27848
0.28446
0.29286
Present 3D (M=50, N=3)
0.27491
0.27849
0.28447
0.29287
Present 3D (M=100, N=3)
0.27491
0.27849
0.28447
0.29287
3D [8, 9]
0.27491
0.27849
0.28447
0.29287Table 3
Assessment 2 (part II). Simply supported aluminum cylinder with Poisson ratioν=0.3 and thickness ratio h/Rα=0.18. Lowest exact natural circular frequencies in no-dimensional form ω-=ω(h/π)ρ/G for imposed axial half-wave number n=1 and several half-wave numbers m in circumferential direction. Comparison between present three-dimensional analysis and three-dimensional analysis by Armenakas et al. [8] as adapted in Bhimaraddi [9]. M is the number of fictitious layers and N is the order of expansion for the exponential matrix.
h
/
R
α
=
0.18
m
=
2
m
=
4
m
=
6
m
=
8
b
/
R
α
=
2
Present 3D (M=1, N=3)
0.05613
0.03607
0.03241
0.02570
Present 3D (M=10, N=3)
0.05652
0.03930
0.04998
0.07822
Present 3D (M=50, N=3)
0.05652
0.03929
0.04996
0.07821
Present 3D (M=100, N=3)
0.05652
0.03929
0.04996
0.07821
3D [8, 9]
0.05652
0.03929
0.04996
0.07821
b
/
R
α
=
1
Present 3D (M=1, N=3)
0.08833
0.07237
0.05717
0.01965
Present 3D (M=10, N=3)
0.09403
0.08546
0.09094
0.11205
Present 3D (M=50, N=3)
0.09402
0.08545
0.09093
0.11205
Present 3D (M=100, N=3)
0.09402
0.08545
0.09093
0.11205
3D [8, 9]
0.09402
0.08545
0.09093
0.11205
b
/
R
α
=
0.5
Present 3D (M=1, N=3)
0.36506
0.37946
0.40165
0.43001
Present 3D (M=10, N=3)
0.18895
0.19467
0.20616
0.22449
Present 3D (M=50, N=3)
0.18894
0.19467
0.20616
0.22450
Present 3D (M=100, N=3)
0.18894
0.19467
0.20616
0.22450
3D [8, 9]
0.18894
0.19467
0.20616
0.22450
b
/
R
α
=
0.25
Present 3D (M=1, N=3)
0.72230
0.72911
0.74032
0.75570
Present 3D (M=10, N=3)
0.50335
0.50934
0.51931
0.53321
Present 3D (M=50, N=3)
0.50338
0.50937
0.51934
0.53325
Present 3D (M=100, N=3)
0.50338
0.50937
0.51934
0.53325
3D [8, 9]
0.50338
0.50937
0.51934
0.53325The assessment proposed by Soldatos and Ye [7] gives three-dimensional results for multilayered orthotropic cylinders(-θ/+θ)j. 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 (E1/E2=E1/E3=40, G12/E2=G13/E2=0.6, G23/E2=0.5, ν12=ν13=ν23=0.25, and homogeneous mass density ρ). The length b is 10 and the radius of curvature Rα equals 10. The circular dimension is a=2πRα. The cylinder is very thick (h=3); that means thickness ratio is h/Rα=0.3. In Table 4 the first three circular frequencies in no-dimensional form ω-=ω(R/π)ρ/G12 are given for half-wave numbers m=0 and n=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 M=100 fictitious layers and N=3 order of expansion. (M=100, N=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.Table 4
Assessment 3. Simply supported orthotropic cylinder with thickness ratioh/Rα=0.3 and orthotropic angle θ equal to 0° or 90°. First three exact natural circular frequencies in no-dimensional form ω-=ω(R/π)ρ/G12 for imposed half-wave numbers m=0 and n=1. Comparison between present three-dimensional analysis and three-dimensional analysis by Soldatos and Ye [7]. M is the number of fictitious layers and N is the order of expansion for the exponential matrix.
Mode
I
II
III
I
II
III
(m=0,  n=1), θ=0°
(m=0, n=1), θ=90°
Present 3D (M=50, N=3)
1.02407
2.55671
4.66395
0.20406
0.49393
0.91564
Present 3D (M=60, N=3)
0.32274
0.41827
3.29830
0.30493
0.81346
7.51166
Present 3D (M=80, N=3)
0.99999
1.29840
1.85418
0.95308
1.00000
3.04558
Present 3D (M=90, N=3)
1.00000
1.29714
2.48260
0.95308
1.00000
3.25249
Present 3D (M=100, N=3)
1.00000
1.29714
2.48260
0.95308
1.00000
3.25249
3D [7]
1.00000
1.29730
2.48260
0.95307
1.00000
3.25250All the new benchmarks proposed in Section4.2 will use (M=1, N=18) solution for plate cases and (M=100, N=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 Figure3. All these structures are simply supported. The square plate has dimensions a=b=1 and thickness ratios a/h=100,50,10,5. The rectangular plate has dimensions a=1 and b=3a and the thickness ratios are a/h=100,50,10,5. The cylindrical shell panel has a radius of curvature Rα=10 and an infinite radius of curvature Rβ in β direction. The dimensions are a=(π/3)Rα and b=20, and the thickness ratios are Rα/h=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 a=2πRα. The other dimension is b=100 and the thickness ratios are Rα/h=1000,100,10,5. The last geometry is the spherical shell panel with radii of curvature Rα=Rβ=10, dimensions a=b=(π/3)Rα, and thickness ratios Rα/h=1000,100,10,5. Each geometry includes two materials: one isotropic layer in aluminium alloy Al2024 with Young modulus E=73GPa, Poisson ratio ν=0.3, and mass density ρ=2800kg/m3 and one orthotropic layer in Gr/Ep with Young modulus components E1=132.38GPa and E2=E3=10.756GPa, shear modulus components G12=G13=5.6537GPa and G23=3.603GPa, Poisson ratio components ν12=ν13=0.24 and ν23=0.49, and mass density ρ=1600kg/m3. The benchmarks are one-layered isotropic square plate (see 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 ω-=ω(a2/h)ρ/E2 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 u*=u/|umax⁡|, v*=v/|vmax⁡|, w*=w/|wmax⁡|, and z*=z~/h.Table 5
Benchmark 1. Simply supported isotropic square plate. First three vibration modes in term of no-dimensional circular frequencyω-=ω(a2/h)ρ/E for several combinations of half-wave numbers m and n. M=1 fictitious layers and N=18 order of expansion for the exponential matrix.
a
/
h
100
50
10
5
(m=0,  n=1)
I mode
2.9861
2.9846
2.9360
2.8012
II mode
194.83
97.417
19.483
9.7417
III mode
329.33
164.66
32.908
16.415
(m=1, n=0)
I mode
2.9861
2.9846
2.9360
2.8012
II mode
194.83
97.417
19.483
9.7417
III mode
329.33
164.66
32.908
16.415
(m=1, n=1)
I mode
5.9713
5.9650
5.7769
5.3036
II mode
275.54
137.77
27.554
13.777
III mode
465.73
232.86
46.502
23.136
(m=1, n=2)
I mode
14.920
14.881
13.805
11.645
II mode
435.66
217.83
43.566
21.783
III mode
736.37
368.14
73.350
36.157
(m=2, n=1)
I mode
14.920
14.881
13.805
11.645
II mode
435.66
217.83
43.566
21.783
III mode
736.37
368.14
73.350
36.157
(m=2, n=2)
I mode
23.860
23.761
21.214
16.882
II mode
551.07
275.54
55.107
27.554
III mode
931.42
465.63
92.545
45.071Table 6
Benchmark 2. Simply supported orthotropic square plate. First three vibration modes in term of no-dimensional circular frequencyω-=ω(a2/h)ρ/E2 for several combinations of half-wave numbers m and n. M=1 fictitious layers and N=18 order of expansion for the exponential matrix.
a
/
h
100
50
10
5
100
50
10
5
(
m
=
0,  n=1),θ=0°
(
m
=
0
,
n
=
1
),  θ=90°
I mode
2.8553
2.8537
2.8058
2.6731
10.007
9.9712
9.0075
7.2059
II mode
227.77
113.88
22.777
11.388
227.77
113.88
22.777
11.388
III mode
314.89
157.44
31.457
15.679
1104.7
552.35
110.41
46.861
(
m
=
1,  n=0),θ=0°
(
m
=
1
,
n
=
0
)
,
θ
=
9
0
°
I mode
10.007
9.9712
9.0075
7.2059
2.8553
2.8537
2.8058
2.6731
II mode
227.77
113.88
22.777
11.388
227.77
113.88
22.777
11.388
III mode
1104.7
552.35
110.41
46.861
314.89
157.44
31.457
15.679
(
m
=
1
,
n
=
1
)
,
θ
=
0
°
(
m
=
1
,
n
=
1
)
,
θ
=
9
0
°
I mode
11.367
11.325
10.226
8.2183
11.367
11.325
10.226
8.2183
II mode
382.03
191.01
38.169
19.032
382.03
191.01
38.169
19.032
III mode
1130.2
565.09
112.94
49.636
1130.2
565.09
112.94
49.636
(
m
=
1
,
n
=
2
)
,
θ
=
0
°
(
m
=
1
,
n
=
2
)
,
θ
=
9
0
°
I mode
17.719
17.648
15.817
12.656
41.011
40.436
29.614
19.220
II mode
652.40
326.17
65.004
32.100
549.29
274.64
54.882
27.380
III mode
1204.5
602.22
120.31
56.667
2222.2
1111.0
190.33
53.344
(
m
=
2
,
n
=
1
)
,
θ
=
0
°
(
m
=
2
,
n
=
1
)
,
θ
=
9
0
∘
I mode
41.011
40.436
29.614
19.220
17.719
17.648
15.817
12.656
II mode
549.29
274.64
54.882
27.380
652.40
326.17
65.004
32.100
III mode
2222.2
1111.0
190.33
53.344
1204.5
602.22
120.31
56.667
(
m
=
2
,
n
=
2
)
,
θ
=
0
°
(
m
=
2
,
n
=
2
)
,
θ
=
9
0
°
I mode
45.302
44.664
32.873
21.768
45.302
44.664
32.873
21.768
II mode
764.05
381.98
76.128
37.625
764.05
381.98
76.128
37.625
III mode
2260.4
1130.1
198.54
59.836
2260.4
1130.1
198.54
59.836Table 7
Benchmark 3. Simply supported isotropic rectangular plate. First three vibration modes in terms of no-dimensional circular frequencyω-=ω(a2/h)ρ/E for several combinations of half-wave numbers m and n. M=1 fictitious layers and N=18 order of expansion for the exponential matrix.
a
/
h
100
50
10
5
(
m
=
0
,
n
=
1
)
I mode
0.3318
0.3318
0.3312
0.3293
II mode
64.944
32.472
6.4944
3.2472
III mode
109.78
54.888
10.977
5.4870
(
m
=
1
,
n
=
0
)
I mode
2.9861
2.9846
2.9360
2.8012
II mode
194.83
97.417
19.483
9.7417
III mode
329.33
164.66
32.908
16.415
(
m
=
1
,
n
=
1
)
I mode
3.3179
3.3159
3.2562
3.0923
II mode
205.37
102.69
20.537
10.269
III mode
347.14
173.56
34.685
17.297
(
m
=
1,  n=2)
I mode
4.3130
4.3097
4.2099
3.9450
II mode
234.16
117.08
23.416
11.708
III mode
395.80
197.89
39.537
19.699
(
m
=
2
,
n
=
1
)
I mode
12.270
12.243
11.497
9.9094
II mode
395.04
197.52
39.504
19.752
III mode
667.72
333.83
66.559
32.908
(
m
=
2, n=2)
I mode
13.264
13.233
12.369
10.572
II mode
410.74
205.37
41.074
20.537
III mode
694.26
347.10
69.186
34.170Table 8
Benchmark 4. Simply supported orthotropic rectangular plate. First three vibration modes in terms of no-dimensional circular frequencyω-=ω(a2/h)ρ/E2 for several combinations of half-wave numbers m and n. M=1 fictitious layers and N=18 order of expansion for the exponential matrix.
a
/
h
100
50
10
5
100
50
10
5
(
m
=
0,  n=1),  θ=0°
(
m
=
0,  n=1),  θ=90°
I mode
0.3173
0.3173
0.3167
0.3148
1.1130
1.1126
1.0987
1.0586
II mode
75.922
37.961
7.5922
3.7961
75.922
37.961
7.5922
3.7961
III mode
104.98
52.483
10.495
5.2459
368.24
184.12
36.822
18.408
(
m
=
1,  n=0),  θ=0°
(
m
=
1,  n=0),  θ=90°
I mode
10.007
9.9712
9.0075
7.2059
2.8553
2.8537
2.8058
2.6731
II mode
227.77
113.88
22.777
11.388
227.77
113.68
22.777
11.388
III mode
1104.7
552.35
110.41
46.861
314.89
157.44
31.457
15.679
(
m
=
1,  n=1), θ=0°
(
m
=
1,  n=1),  θ=90°
I mode
10.127
10.091
9.1126
7.2879
3.4236
3.4213
3.3511
3.1611
II mode
249.70
124.85
24.967
12.480
312.77
156.38
31.251
15.586
III mode
1107.6
553.78
110.69
47.184
441.11
220.55
44.099
22.031
(
m
=
1,  n=2),  θ=0°
(
m
=
1,  n=2),  θ=90°
I mode
10.538
10.500
9.4774
7.5865
6.1047
6.0947
5.8032
5.1319
II mode
306.16
153.08
30.604
15.284
341.88
170.94
34.157
17.031
III mode
1116.1
558.04
111.54
48.130
774.35
387.17
77.407
38.652
(
m
=
2,  n=1),  θ=0°
(
m
=
2,  n=1),  θ=90°
I mode
40.001
39.445
28.902
18.694
11.868
11.841
11.093
9.5179
II mode
466.89
233.44
46.684
23.337
560.57
280.28
56.021
27.943
III mode
2210.8
1105.3
187.77
51.133
656.59
328.31
65.435
32.318
(
m
=
2,  n=2),  θ=0°
(
m
=
2,  n=2),  θ=90°
I mode
40.364
39.801
29.152
18.872
13.685
13.649
12.644
10.634
II mode
499.39
249.69
49.921
24.938
625.52
312.73
62.343
30.811
III mode
2215.1
1107.5
188.74
51.982
882.21
441.09
88.126
43.879Table 9
Benchmark 5. Simply supported isotropic cylindrical shell panel. First three vibration modes in terms of no-dimensional circular frequencyω-=ω(a2/h)ρ/E for several combinations of half-wave numbers m and n. M=100 fictitious layers and N=3 order of expansion for the exponential matrix.
R
α
/
h
1000
100
10
5
(
m
=
0
,
n
=
1
)
I mode
1067.1
106.71
10.683
5.3415
II mode
1068.3
106.83
10.690
5.3727
III mode
1855.7
185.57
18.557
9.2797
(
m
=
1
,
n
=
0
)
I mode
2.5478
2.5184
2.4853
2.3948
II mode
2040.3
204.03
20.411
10.216
III mode
3635.3
363.51
36.229
17.917
(
m
=
1
,
n
=
1
)
I mode
222.72
22.503
3.8758
3.2156
II mode
2364.5
236.45
23.656
11.842
III mode
4017.4
401.73
40.060
19.846
(
m
=
1
,
n
=
2
)
I mode
552.89
55.574
7.7851
5.7590
II mode
3031.5
303.16
30.328
15.181
III mode
5048.7
504.86
50.367
24.983
(
m
=
2
,
n
=
1
)
I mode
70.382
14.052
11.533
10.021
II mode
4230.8
423.08
42.325
21.174
III mode
7214.2
721.38
71.724
35.120
(
m
=
2
,
n
=
2
)
I mode
232.90
27.451
13.793
11.643
II mode
4639.2
463.92
46.414
23.226
III mode
7847.2
784.68
78.020
38.182Table 10
Benchmark 6. Simply supported orthotropic cylindrical shell panel. First three vibration modes in term of no-dimensional circular frequencyω-=ω(a2/h)ρ/E2 for several combinations of half-wave numbers m and n. M=100 fictitious layers and N=3 order of expansion for the exponential matrix.
R
α
/
h
1000
100
10
5
1000
100
10
5
(
m
=
0
,
n
=
1
)
,
θ
=
0
°
(
m
=
0
,
n
=
1
)
,
θ
=
9
0
°
I mode
1248.9
124.89
12.489
6.2444
1096.5
109.67
11.281
6.0108
II mode
1721.6
172.16
17.230
8.6360
1248.9
124.89
12.489
6.2444
III mode
3858.5
385.85
38.590
19.293
6057.8
605.78
60.570
30.272
(
m
=
1
,
n
=
0
)
,
θ
=
0
°
(
m
=
1
,
n
=
0
)
,
θ
=
9
0
°
I mode
8.4450
8.4392
7.6484
6.1872
2.3272
2.4079
2.3754
2.2863
II mode
2385.2
238.52
23.861
11.941
2385.2
238.52
23.861
11.943
III mode
12194
1219.4
121.48
51.358
3476.0
347.58
34.633
17.114
(
m
=
1
,
n
=
1
)
,
θ
=
0
°
(
m
=
1
,
n
=
1
)
,
θ
=
9
0
°
I mode
230.37
24.568
8.0953
6.3709
352.48
35.490
5.3604
4.1227
II mode
2939.9
294.00
29.454
14.782
3588.6
358.85
35.788
17.735
III mode
12257
1225.7
122.10
52.106
6559.0
655.90
65.574
32.756
(
m
=
1
,
n
=
2
)
,
θ
=
0
°
(
m
=
1
,
n
=
2
)
,
θ
=
9
0
°
I mode
641.70
64.887
10.877
6.8670
644.27
65.496
12.428
9.0851
II mode
4155.7
415.58
41.604
20.824
4164.3
416.43
41.609
20.709
III mode
12445
1244.4
123.94
54.240
12370
1237.0
123.62
54.652
(
m
=
2
,
n
=
1
)
,
θ
=
0
°
(
m
=
2
,
n
=
1
)
,
θ
=
9
0
°
I mode
79.744
39.112
28.352
18.587
160.90
20.150
11.545
9.9019
II mode
5062.8
506.29
50.681
25.371
6476.9
647.66
64.461
31.673
III mode
23492
2349.1
205.91
56.149
7983.5
798.33
79.684
39.571
(
m
=
2
,
n
=
2
)
,
θ
=
0
°
(
m
=
2
,
n
=
2
)
,
θ
=
9
0
°
I mode
243.62
46.084
29.034
19.077
364.43
40.611
16.527
13.080
II mode
5849.5
584.96
58.566
29.293
6951.1
695.08
69.132
33.884
III mode
23598
2359.6
208.24
58.169
13115
1311.5
131.01
62.003Table 11
Benchmark 7. Simply supported isotropic cylinder. First three vibration modes in term of no-dimensional circular frequencyω-=ω(a2/h)ρ/E for several combinations of half-wave numbers m and n. M=100 fictitious layers and N=3 order of expansion for the exponential matrix.
R
α
/
h
1000
100
10
5
(m=0, n=1)
I mode
7691.7
769.17
76.917
38.459
II mode
12342
1234.2
123.42
61.710
III mode
41586
4158.6
416.30
208.81
(
m
=
2, n=1)
I mode
2445.9
244.59
24.492
12.295
II mode
27270
2727.0
272.81
136.56
III mode
59189
5918.8
591.08
294.28
(
m
=
2,  n=2)
I mode
7712.1
771.23
77.252
38.820
II mode
33375
3337.5
333.85
167.09
III mode
61353
6135.3
612.74
305.15
(
m
=
2,  n=3)
I mode
13536
1353.7
135.69
68.323
II mode
40081
4008.1
400.89
200.56
III mode
65490
6548.9
654.17
325.96Table 12
Benchmark 8. Simply supported orthotropic cylinder. First three vibration modes in terms of no-dimensional circular frequencyω-=ω(a2/h)ρ/E2 for several combinations of half-wave numbers m and n. M=100 fictitious layers and N=3 order of expansion for the exponential matrix.
R
α
/
h
1000
100
10
5
1000
100
10
5
(
m
=
0,  n=1), θ=0°
(
m
=
0, n=1), θ=90°
I mode
8991.9
899.19
89.919
44.960
8991.9
899.19
89.919
44.959
II mode
12402
1240.2
124.16
62.281
39097
3909.7
391.39
196.33
III mode
138826
13883
1389.6
696.66
44039
4403.9
440.45
220.32
(
m
=
2, n=1), θ=0°
(
m
=
2, n=1), θ=90°
I mode
2486.4
248.65
24.921
12.541
5295.1
529.53
53.084
26.738
II mode
31697
3169.8
317.14
158.81
50390
5039.0
503.83
251.80
III mode
196438
19643
1961.3
974.72
58023
5802.3
579.62
288.86
(
m
=
2,  n=2),  θ=0°
(
m
=
2, n=2),  θ=90°
I mode
8071.9
807.21
80.903
40.715
11839
1184.0
119.10
60.562
II mode
38989
3898.9
390.16
195.44
56607
5660.7
565.47
281.81
III mode
196773
19677
1964.6
976.14
92402
9240.2
924.00
461.97
(
m
=
2, n=3),  θ=0°
(
m
=
2,  n=3), θ=90°
I mode
14640
1464.1
146.76
73.903
17656
1765.9
178.90
92.578
II mode
48229
4823.0
482.66
241.82
58769
5876.8
587.41
293.24
III mode
197336
19733
1970.0
978.55
134283
13428
1342.8
671.29Table 13
Benchmark 9. Simply supported isotropic spherical shell panel. First three vibration modes in terms of no-dimensional circular frequencyω-=ω(a2/h)ρ/E for several combinations of half-wave numbers m and n. M=100 fictitious layers and N=3 order of expansion for the exponential matrix.
R
α
/
h
1000
100
10
5
(
m
=
0, n=1)
I mode
999.01
99.924
10.213
5.3986
II mode
2040.3
204.02
20.360
10.115
III mode
3785.7
378.55
37.677
18.556
(
m
=
1,  n=0)
I mode
999.01
99.924
10.213
5.3986
II mode
2040.3
204.02
20.360
10.116
III mode
3785.7
378.55
37.677
18.556
(
m
=
1,  n=1)
I mode
1046.3
104.76
11.564
6.9554
II mode
2885.4
288.53
28.793
14.302
III mode
5111.6
511.13
50.821
24.937
(
m
=
1, n=2)
I mode
1076.3
108.52
16.906
12.400
II mode
4562.2
456.21
45.524
22.598
III mode
7857.9
785.72
77.933
37.847
(
m
=
2, n=1)
I mode
1076.3
108.52
16.906
12.400
II mode
4562.2
456.21
45.524
22.598
III mode
7857.9
785.72
77.933
37.847
(
m
=
2,  n=2)
I mode
1084.1
110.77
23.273
17.513
II mode
5770.8
577.07
57.581
28.565
III mode
9869.7
986.87
97.645
46.869Table 14
Benchmark 10. Simply supported orthotropic spherical shell panel. First three vibration modes in terms of no-dimensional circular frequencyω-=ω(a2/h)ρ/E2 for several combinations of half-wave numbers m and n. M=100 fictitious layers and N=3 order of expansion for the exponential matrix.
R
α
/
h
1000
100
10
5
1000
100
10
5
(
m
=
0, n=1), θ=0∘
(
m
=
0,  n=1), θ=90∘
I mode
2385.2
238.51
23.802
11.825
1037.9
104.08
12.585
7.7554
II mode
2902.0
290.20
29.026
14.513
2385.2
238.51
23.802
11.824
III mode
4371.6
437.15
43.559
21.518
12223
1222.3
121.65
51.888
(
m
=
1, n=0), θ=0∘
(
m
=
1, n=0), θ=90∘
I mode
1037.9
104.08
12.583
7.7554
2385.2
238.51
23.802
11.825
II mode
2385.2
238.51
23.802
11.824
2902.0
290.20
29.026
14.513
III mode
12223
1222.3
121.65
51.888
4371.6
437.15
43.559
21.518
(
m
=
1, n=1), θ=0∘
(
m
=
1, n=1), θ=90∘
I mode
1209.7
121.32
14.754
9.1755
1209.7
121.32
14.754
9.1755
II mode
4016.5
401.63
40.041
19.819
4016.5
401.63
40.041
19.819
III mode
12456
1245.5
123.88
54.506
12456
1245.5
123.88
54.506
(
m
=
1,  n=2), θ=0∘
(
m
=
1, n=2), θ=90∘
I mode
1682.9
169.01
21.986
14.305
1114.0
117.90
30.685
19.645
II mode
6838.7
683.84
68.035
33.373
5758.3
575.82
57.401
28.398
III mode
13147
1314.6
130.57
60.847
23599
2359.8
208.50
58.572
(
m
=
2,  n=1), θ=0∘
(
m
=
2, n=1), θ=90∘
I mode
1114.0
117.90
30.685
19.645
1682.9
169.01
21.986
14.305
II mode
5758.3
575.82
57.401
28.398
6838.7
683.84
68.035
33.373
III mode
23599
2359.8
208.50
58.572
13147
1314.6
130.57
60.847
(
m
=
2, n=2), θ=0∘
(
m
=
2, n=2), θ=90∘
I mode
1260.9
133.19
34.372
22.411
1260.9
133.19
34.372
22.411
II mode
8009.4
800.89
79.625
39.004
8009.4
800.89
79.625
39.004
III mode
23985
2398.3
216.62
65.349
23985
2398.3
216.62
65.349Figure 4
Benchmark 1, simply supported isotropic square plate. First three vibration modes in terms of displacement components through the thickness for thickness ratioa/h=10 and several combinations of half-wave numbers.Figure 5
Benchmark 2, simply supported orthotropic square plate(θ=0°). First three vibration modes in terms of displacement components through the thickness for thickness ratio a/h=10 and several combinations of half-wave numbers.Figure 6
Benchmark 2, simply supported orthotropic square plate (θ=90°). First three vibration modes in terms of displacement components through the thickness for thickness ratio a/h=10 and several combinations of half-wave numbers.Figure 7
Benchmark 7, simply supported isotropic cylinder. First three vibration modes in terms of displacement components through the thickness for thickness ratioRα/h=10 and several combinations of half-wave numbers.Figure 8
Benchmark 8, simply supported orthotropic cylinder(θ=0°). First three vibration modes in terms of displacement components through the thickness for thickness ratio Rα/h=10 and several combinations of half-wave numbers.Figure 9
Benchmark 8, simply supported orthotropic cylinder(θ=90°). First three vibration modes in terms of displacement components through the thickness for thickness ratio Rα/h=10 and several combinations of half-wave numbers.Figure 10
Benchmark 9, simply supported isotropic spherical shell panel. First three vibration modes in terms of displacement components through the thickness for thickness ratioRα/h=5 and several combinations of half-wave numbers.Figure 11
Benchmark 10, simply supported orthotropic spherical shell panel(θ=0°). First three vibration modes in terms of displacement components through the thickness for thickness ratio Rα/h=5 and several combinations of half-wave numbers.Figure 12
Benchmark 10, simply supported orthotropic spherical shell panel(θ=90°). First three vibration modes in terms of displacement components through the thickness for thickness ratio Rα/h=5 and several combinations of half-wave numbers.Table5 presents a square isotropic plate (benchmark 1). It shows the first three vibration modes for each possible combination of half-wave numbers (m,n) 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 (a/h=10) in terms of the three displacement components u*, v*, and w* 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 w* through the thickness is zero. The first mode has a constant transverse displacement w* through the thickness and linear in-plane displacements u* 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 n=0. The third vibration mode always has linear transverse displacement w* through the thickness for each pair of half-wave numbers. In-plane displacements u* 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 n.Table6 shows a square orthotropic plate (benchmark 2) and lists the first three vibration modes for each possible combination of half-wave numbers (m,n) 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 u*, v*, and w* for half-wave numbers (1,0), (1,1), (2,1), and (2,2) and for a thick plate (a/h=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 half-wave 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.Tables7 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 a and b are not the same (b=3a). 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).Table9 gives the first three vibration modes of the isotropic cylindrical shell panel (benchmark 5) for several combinations of half-wave numbers (m,n). The shell is isotropic but the frequencies are always different for each couple of (m,n) because of the geometry with different a and b dimensions. The same considerations are also valid for benchmark 6 in Table 10 (the orthotropic cylindrical shell panel). The geometry with different a and b 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.Tables11 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 (m,n) always has constant displacements u* and w*, and a quasi-constant v* displacement, but in the case m=0 and n=1 the only component different from zero is u*. The second mode always has constant displacements v* and w* and a linear u* displacement, but in the case m=0 and n=1u* is zero and w* and v* are constant. The third mode has quasi-constant u*, v*, and w* displacements, but in the case m=0 and n=1 the component u* 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 (m,n) 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 w* 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∘).Tables13 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 a and b. 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 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 theHα effect in cylindrical shell panels and both Hα and Hβ 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 Hα=1 (which means z/Rα≃0) for the benchmark 5 concerning the isotropic cylindrical shell panel. The approximation Hα=1 is valid for very thin shells (Rα/h equals 1000 or 100) even if it could have some problems for particular vibration modes and higher values of half-wave numbers (m,n). For thick shells the use of the exact geometry without the approximation Hα=1 is mandatory. Table 16 shows the comparison between the present 3D solution using the exact geometry and the same solution using the approximation Hα=Hβ=1 (which means z/Rα=z/Rβ≃0) for the benchmark 10 concerning the orthotropic spherical shell panel. The approximation Hα=Hβ=1 is valid for thin shells (Rα/h 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 z/Rα=z/Rβ≃0 is not valid.Table 15
Effects ofHα=(1+z/Rα) terms in benchmark 5 for simply supported isotropic cylindrical shell panel. First three vibration modes in terms of no-dimensional circular frequency ω-=ω(a2/h)ρ/E for several combinations of half-wave numbers m and n. M=100 fictitious layers and N=3 order of expansion for the exponential matrix.
R
α
/
h
1000
100
10
5
(
m
=
0, n=1)
I mode
1067.1
106.71
10.683
5.3415
II mode
1068.3
106.83
10.690
5.3727
III mode
1855.7
185.57
18.557
9.2797
I mode (Hα=1)
1067.1
106.71
10.686
5.3486
II mode (Hα=1)
1068.3
106.83
10.688
5.3684
III mode (Hα=1)
1855.7
185.56
18.550
9.2651
(
m
=
1,  n=0)
I mode
2.5478
2.5184
2.4853
2.3948
II mode
2040.3
204.03
20.411
10.216
III mode
3635.3
363.51
36.229
17.917
I mode (Hα=1)
2.7747
2.7473
2.7029
2.5828
II mode (Hα=1)
2040.3
204.03
20.403
10.201
III mode (Hα=1)
3635.3
363.52
36.303
18.076
(
m
=
1, n=1)
I mode
222.72
22.503
3.8758
3.2156
II mode
2364.5
236.45
23.656
11.842
III mode
4017.4
401.73
40.060
19.846
I mode (Hα=1)
222.73
22.535
4.0440
3.3738
II mode (Hα=1)
2364.5
236.45
23.668
11.868
III mode (Hα=1)
4017.4
401.73
40.115
19.965Table 16
Effects ofHα=(1+z/Rα) and Hβ=(1+z/Rβ) terms in benchmark 10 for simply supported orthotropic spherical shell panel. First three vibration modes in term of no-dimensional circular frequency ω-=ω(a2/h)ρ/E2 for several combinations of half-wave numbers m and n. M=100 fictitious layers and N=3 order of expansion for the exponential matrix.
R
α
/
h
1000
100
10
5
1000
100
10
5
(
m
=
0, n=1), θ=0°
(
m
=
0, n=1), θ=90°
I mode
2385.2
238.51
23.802
11.825
1037.9
104.08
12.585
7.7554
II mode
2902.0
290.20
29.026
14.513
2385.2
238.51
23.802
11.824
III mode
4371.6
437.15
43.559
21.518
12223
1222.3
121.65
51.888
I mode (Hα=Hβ=1)
2385.2
238.52
23.858
11.939
1037.9
104.09
12.602
7.7410
II mode (Hα=Hβ=1)
2902.0
290.20
29.030
14.529
2385.2
238.52
23.856
11.934
III mode (Hα=Hβ=1)
4371.6
437.15
43.623
21.667
12223
1222.3
122.04
51.743
(
m
=
1, n=0), θ=0°
(
m
=
1, n=0), θ=90°
I mode
1037.9
104.08
12.583
7.7554
2385.2
238.51
23.802
11.825
II mode
2385.2
238.51
23.802
11.824
2902.0
290.20
29.026
14.513
III mode
12223
1222.3
121.65
51.888
4371.6
437.15
43.559
21.518
I mode (Hα=Hβ=1)
1037.9
104.09
12.602
7.7410
2385.2
238.52
23.858
11.939
II mode (Hα=Hβ=1)
2385.2
238.52
23.856
11.934
2902.0
290.20
29.030
14.529
III mode (Hα=Hβ=1)
12223
1222.3
122.04
51.743
4371.6
437.15
43.623
21.667
(
m
=
1, n=1), θ=0°
(
m
=
1, n=1), θ=90°
I mode
1209.7
121.32
14.754
9.1755
1209.7
121.32
14.754
9.1755
II mode
4016.5
401.63
40.041
19.819
4016.5
401.63
40.041
19.819
III mode
12456
1245.5
123.88
54.506
12456
1245.5
123.88
54.506
I mode (Hα=Hβ=1)
1209.7
121.33
14.810
9.2169
1209.7
121.33
14.810
9.2169
II mode (Hα=Hβ=1)
4016.5
401.64
40.137
20.022
4016.5
401.64
40.137
20.022
III mode (Hα=Hβ=1)
12456
1245.6
124.29
54.312
12456
1245.6
124.29
54.312


## 4.1. 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 E=70GPa, Poisson ratio ν=0.3, and mass density ρ=2702kg/m3. The thickness value is h=1. Two different thickness ratios are considered, that is, a/h=10 and a/h=10. The first six circular frequencies in no-dimensional form ω-=ω(a2/h)ρ/E are given in Table 1 for imposed half-wave numbers m=n=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 M number of fictitious layers and the N 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 Hα and Hβ. For M=1, an order of expansion N=18 for the exponential matrix always gives the correct solution for each vibration mode and thickness ratio a/h. The calculation of the exponential matrix is not time consuming. An alternative method could be the use of M fictitious layers and a lower order of expansion N=3 for the exponential matrix. In this last case, M=100 fictitious layers always give the correct solution for each vibration mode and thickness ratio. The use of M=100 layers is not time consuming because the E matrix has always 6×6 dimension for any value of M layers, even if the proposed method uses a layerwise approach. Solutions with M=1 and N=18 or with M=100 and N=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.Table 1
Assessment 1. Simply supported aluminum square plate with Poisson ratioν=0.3. First six exact natural circular frequencies in no-dimensional form ω-=ω(a2/h)ρ/E for imposed half-wave numbers m=n=1. Comparison between present three-dimensional analysis and three-dimensional analyses by Vel and Batra [34] and Srinivas et al. [22]. M is the number of fictitious layers and N is the order of expansion for the exponential matrix.
Mode
I
II
III
IV
V
VI
a
/
h
=
10
Present 3D (M=1, N=3)
5.3370
27.554
46.493
154.39
168.09
268.36
Present 3D (M=1, N=7)
5.7769
27.554
46.502
192.91
197.98
329.35
Present 3D (M=1, N=10)
5.7769
27.554
46.502
196.21
201.87
308.34
Present 3D (M=1, N=15)
5.7769
27.554
46.502
196.77
201.34
358.23
Present 3D (M=1, N=18)
5.7769
27.554
46.502
196.77
201.34
357.42
Present 3D (M=10, N=3)
5.7769
27.554
46.502
196.71
201.27
355.40
Present 3D (M=50, N=3)
5.7769
27.554
46.502
196.77
201.34
351.90
Present 3D (M=100, N=3)
5.7769
27.554
46.502
196.77
201.34
357.42
3D [34]
5.7769
27.554
46.503
196.77
201.34
357.42
3D [22]
5.7769
27.554
46.502
196.77
201.34
357.42
a
/
h
=
10
Present 3D (M=1, N=3)
/
8.7132
14.182
17.512
/
/
Present 3D (M=1, N=7)
4.6598
8.7132
14.463
20.987
23.999
29.925
Present 3D (M=1, N=10)
4.6582
8.7132
14.463
21.383
25.009
31.399
Present 3D (M=1, N=15)
4.6581
8.7132
14.463
21.343
24.830
34.049
Present 3D (M=1, N=18)
4.6581
8.7132
14.463
21.343
24.830
33.982
Present 3D (M=10, N=3)
4.6580
8.7132
14.463
21.337
24.821
33.939
Present 3D (M=50, N=3)
4.6581
8.7132
14.463
21.343
24.830
33.923
Present 3D (M=100, N=3)
4.6581
8.7132
14.463
21.343
24.830
33.982
3D [34]
4.6582
8.7132
14.463
21.343
24.830
33.982
3D [22]
4.6582
8.7134
14.462
21.343
24.830
33.982The 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 Rα=1 and infinite radius of curvature Rβ in β direction. The circular dimension is a=2πRα; the thickness h is 0.12 and 0.18; that means thickness ratios are h/Rα=0.12 and h/Rα=0.18, respectively. Several ratios cylinder length/radius of curvature (b/Rα) are investigated in the case of b=2,1,0.5,0.25. The fundamental circular frequency in no-dimensional form ω-=ω(h/π)ρ/G is given in Tables 2 and 3 (thickness ratio h/Rα equals 0.12 ad 0.18, resp.) for axial half-wave number n=1 and several circular half-wave numbers m=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 Hα and Hβ. The convergence study is shown for a value N=3 for the expansion of the exponential matrix. M=100 fictitious layers always give the correct solution for each thickness ratio h/Rα (see Tables 2 and 3), length b of the cylinder, and half-wave numbers m in circumferential direction. The present method is not time consuming; M=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.Table 2
Assessment 2 (part I). Simply supported aluminum cylinder with Poisson ratioν=0.3 and thickness ratio h/Rα=0.12. Lowest exact natural circular frequencies in no-dimensional form ω-=ω(h/π)ρ/G for imposed axial half-wave number n=1 and several half-wave numbers m in circumferential direction. Comparison between present three-dimensional analysis and three-dimensional analysis by Armenakas et al. [8] as adapted in Bhimaraddi [9]. M is the number of fictitious layers and N is the order of expansion for the exponential matrix.
h
/
R
α
=
0.12
m
=
2
m
=
4
m
=
6
m
=
8
b
/
R
α
=
2
Present 3D (M=1, N=3)
0.03724
0.02323
0.02175
0.02846
Present 3D (M=10, N=3)
0.03730
0.02360
0.02464
0.03687
Present 3D (M=50, N=3)
0.03730
0.02359
0.02463
0.03686
Present 3D (M=100, N=3)
0.03730
0.02359
0.02462
0.03686
3D [8, 9]
0.03730
0.02359
0.02462
0.03686
b
/
R
α
=
1
Present 3D (M=1, N=3)
0.05769
0.04772
0.04229
0.04272
Present 3D (M=10, N=3)
0.05853
0.04979
0.04790
0.05546
Present 3D (M=50, N=3)
0.05853
0.04978
0.04789
0.05545
Present 3D (M=100, N=3)
0.05853
0.04978
0.04789
0.05545
3D [8, 9]
0.05853
0.04978
0.04789
0.05545
b
/
R
α
=
0.5
Present 3D (M=1, N=3)
0.08017
0.07732
0.07308
0.06667
Present 3D (M=10, N=3)
0.10057
0.10234
0.10688
0.11528
Present 3D (M=50, N=3)
0.10057
0.10234
0.10688
0.11528
Present 3D (M=100, N=3)
0.10057
0.10234
0.10688
0.11528
3D [8, 9]
0.10057
0.10234
0.10688
0.11528
b
/
R
α
=
0.25
Present 3D (M=1, N=3)
0.48156
0.48619
0.49378
0.50419
Present 3D (M=10, N=3)
0.27489
0.27848
0.28446
0.29286
Present 3D (M=50, N=3)
0.27491
0.27849
0.28447
0.29287
Present 3D (M=100, N=3)
0.27491
0.27849
0.28447
0.29287
3D [8, 9]
0.27491
0.27849
0.28447
0.29287Table 3
Assessment 2 (part II). Simply supported aluminum cylinder with Poisson ratioν=0.3 and thickness ratio h/Rα=0.18. Lowest exact natural circular frequencies in no-dimensional form ω-=ω(h/π)ρ/G for imposed axial half-wave number n=1 and several half-wave numbers m in circumferential direction. Comparison between present three-dimensional analysis and three-dimensional analysis by Armenakas et al. [8] as adapted in Bhimaraddi [9]. M is the number of fictitious layers and N is the order of expansion for the exponential matrix.
h
/
R
α
=
0.18
m
=
2
m
=
4
m
=
6
m
=
8
b
/
R
α
=
2
Present 3D (M=1, N=3)
0.05613
0.03607
0.03241
0.02570
Present 3D (M=10, N=3)
0.05652
0.03930
0.04998
0.07822
Present 3D (M=50, N=3)
0.05652
0.03929
0.04996
0.07821
Present 3D (M=100, N=3)
0.05652
0.03929
0.04996
0.07821
3D [8, 9]
0.05652
0.03929
0.04996
0.07821
b
/
R
α
=
1
Present 3D (M=1, N=3)
0.08833
0.07237
0.05717
0.01965
Present 3D (M=10, N=3)
0.09403
0.08546
0.09094
0.11205
Present 3D (M=50, N=3)
0.09402
0.08545
0.09093
0.11205
Present 3D (M=100, N=3)
0.09402
0.08545
0.09093
0.11205
3D [8, 9]
0.09402
0.08545
0.09093
0.11205
b
/
R
α
=
0.5
Present 3D (M=1, N=3)
0.36506
0.37946
0.40165
0.43001
Present 3D (M=10, N=3)
0.18895
0.19467
0.20616
0.22449
Present 3D (M=50, N=3)
0.18894
0.19467
0.20616
0.22450
Present 3D (M=100, N=3)
0.18894
0.19467
0.20616
0.22450
3D [8, 9]
0.18894
0.19467
0.20616
0.22450
b
/
R
α
=
0.25
Present 3D (M=1, N=3)
0.72230
0.72911
0.74032
0.75570
Present 3D (M=10, N=3)
0.50335
0.50934
0.51931
0.53321
Present 3D (M=50, N=3)
0.50338
0.50937
0.51934
0.53325
Present 3D (M=100, N=3)
0.50338
0.50937
0.51934
0.53325
3D [8, 9]
0.50338
0.50937
0.51934
0.53325The assessment proposed by Soldatos and Ye [7] gives three-dimensional results for multilayered orthotropic cylinders(-θ/+θ)j. 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 (E1/E2=E1/E3=40, G12/E2=G13/E2=0.6, G23/E2=0.5, ν12=ν13=ν23=0.25, and homogeneous mass density ρ). The length b is 10 and the radius of curvature Rα equals 10. The circular dimension is a=2πRα. The cylinder is very thick (h=3); that means thickness ratio is h/Rα=0.3. In Table 4 the first three circular frequencies in no-dimensional form ω-=ω(R/π)ρ/G12 are given for half-wave numbers m=0 and n=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 M=100 fictitious layers and N=3 order of expansion. (M=100, N=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.Table 4
Assessment 3. Simply supported orthotropic cylinder with thickness ratioh/Rα=0.3 and orthotropic angle θ equal to 0° or 90°. First three exact natural circular frequencies in no-dimensional form ω-=ω(R/π)ρ/G12 for imposed half-wave numbers m=0 and n=1. Comparison between present three-dimensional analysis and three-dimensional analysis by Soldatos and Ye [7]. M is the number of fictitious layers and N is the order of expansion for the exponential matrix.
Mode
I
II
III
I
II
III
(m=0,  n=1), θ=0°
(m=0, n=1), θ=90°
Present 3D (M=50, N=3)
1.02407
2.55671
4.66395
0.20406
0.49393
0.91564
Present 3D (M=60, N=3)
0.32274
0.41827
3.29830
0.30493
0.81346
7.51166
Present 3D (M=80, N=3)
0.99999
1.29840
1.85418
0.95308
1.00000
3.04558
Present 3D (M=90, N=3)
1.00000
1.29714
2.48260
0.95308
1.00000
3.25249
Present 3D (M=100, N=3)
1.00000
1.29714
2.48260
0.95308
1.00000
3.25249
3D [7]
1.00000
1.29730
2.48260
0.95307
1.00000
3.25250All the new benchmarks proposed in Section4.2 will use (M=1, N=18) solution for plate cases and (M=100, N=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 Figure3. All these structures are simply supported. The square plate has dimensions a=b=1 and thickness ratios a/h=100,50,10,5. The rectangular plate has dimensions a=1 and b=3a and the thickness ratios are a/h=100,50,10,5. The cylindrical shell panel has a radius of curvature Rα=10 and an infinite radius of curvature Rβ in β direction. The dimensions are a=(π/3)Rα and b=20, and the thickness ratios are Rα/h=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 a=2πRα. The other dimension is b=100 and the thickness ratios are Rα/h=1000,100,10,5. The last geometry is the spherical shell panel with radii of curvature Rα=Rβ=10, dimensions a=b=(π/3)Rα, and thickness ratios Rα/h=1000,100,10,5. Each geometry includes two materials: one isotropic layer in aluminium alloy Al2024 with Young modulus E=73GPa, Poisson ratio ν=0.3, and mass density ρ=2800kg/m3 and one orthotropic layer in Gr/Ep with Young modulus components E1=132.38GPa and E2=E3=10.756GPa, shear modulus components G12=G13=5.6537GPa and G23=3.603GPa, Poisson ratio components ν12=ν13=0.24 and ν23=0.49, and mass density ρ=1600kg/m3. The benchmarks are one-layered isotropic square plate (see 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 ω-=ω(a2/h)ρ/E2 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 u*=u/|umax⁡|, v*=v/|vmax⁡|, w*=w/|wmax⁡|, and z*=z~/h.Table 5
Benchmark 1. Simply supported isotropic square plate. First three vibration modes in term of no-dimensional circular frequencyω-=ω(a2/h)ρ/E for several combinations of half-wave numbers m and n. M=1 fictitious layers and N=18 order of expansion for the exponential matrix.
a
/
h
100
50
10
5
(m=0,  n=1)
I mode
2.9861
2.9846
2.9360
2.8012
II mode
194.83
97.417
19.483
9.7417
III mode
329.33
164.66
32.908
16.415
(m=1, n=0)
I mode
2.9861
2.9846
2.9360
2.8012
II mode
194.83
97.417
19.483
9.7417
III mode
329.33
164.66
32.908
16.415
(m=1, n=1)
I mode
5.9713
5.9650
5.7769
5.3036
II mode
275.54
137.77
27.554
13.777
III mode
465.73
232.86
46.502
23.136
(m=1, n=2)
I mode
14.920
14.881
13.805
11.645
II mode
435.66
217.83
43.566
21.783
III mode
736.37
368.14
73.350
36.157
(m=2, n=1)
I mode
14.920
14.881
13.805
11.645
II mode
435.66
217.83
43.566
21.783
III mode
736.37
368.14
73.350
36.157
(m=2, n=2)
I mode
23.860
23.761
21.214
16.882
II mode
551.07
275.54
55.107
27.554
III mode
931.42
465.63
92.545
45.071Table 6
Benchmark 2. Simply supported orthotropic square plate. First three vibration modes in term of no-dimensional circular frequencyω-=ω(a2/h)ρ/E2 for several combinations of half-wave numbers m and n. M=1 fictitious layers and N=18 order of expansion for the exponential matrix.
a
/
h
100
50
10
5
100
50
10
5
(
m
=
0,  n=1),θ=0°
(
m
=
0
,
n
=
1
),  θ=90°
I mode
2.8553
2.8537
2.8058
2.6731
10.007
9.9712
9.0075
7.2059
II mode
227.77
113.88
22.777
11.388
227.77
113.88
22.777
11.388
III mode
314.89
157.44
31.457
15.679
1104.7
552.35
110.41
46.861
(
m
=
1,  n=0),θ=0°
(
m
=
1
,
n
=
0
)
,
θ
=
9
0
°
I mode
10.007
9.9712
9.0075
7.2059
2.8553
2.8537
2.8058
2.6731
II mode
227.77
113.88
22.777
11.388
227.77
113.88
22.777
11.388
III mode
1104.7
552.35
110.41
46.861
314.89
157.44
31.457
15.679
(
m
=
1
,
n
=
1
)
,
θ
=
0
°
(
m
=
1
,
n
=
1
)
,
θ
=
9
0
°
I mode
11.367
11.325
10.226
8.2183
11.367
11.325
10.226
8.2183
II mode
382.03
191.01
38.169
19.032
382.03
191.01
38.169
19.032
III mode
1130.2
565.09
112.94
49.636
1130.2
565.09
112.94
49.636
(
m
=
1
,
n
=
2
)
,
θ
=
0
°
(
m
=
1
,
n
=
2
)
,
θ
=
9
0
°
I mode
17.719
17.648
15.817
12.656
41.011
40.436
29.614
19.220
II mode
652.40
326.17
65.004
32.100
549.29
274.64
54.882
27.380
III mode
1204.5
602.22
120.31
56.667
2222.2
1111.0
190.33
53.344
(
m
=
2
,
n
=
1
)
,
θ
=
0
°
(
m
=
2
,
n
=
1
)
,
θ
=
9
0
∘
I mode
41.011
40.436
29.614
19.220
17.719
17.648
15.817
12.656
II mode
549.29
274.64
54.882
27.380
652.40
326.17
65.004
32.100
III mode
2222.2
1111.0
190.33
53.344
1204.5
602.22
120.31
56.667
(
m
=
2
,
n
=
2
)
,
θ
=
0
°
(
m
=
2
,
n
=
2
)
,
θ
=
9
0
°
I mode
45.302
44.664
32.873
21.768
45.302
44.664
32.873
21.768
II mode
764.05
381.98
76.128
37.625
764.05
381.98
76.128
37.625
III mode
2260.4
1130.1
198.54
59.836
2260.4
1130.1
198.54
59.836Table 7
Benchmark 3. Simply supported isotropic rectangular plate. First three vibration modes in terms of no-dimensional circular frequencyω-=ω(a2/h)ρ/E for several combinations of half-wave numbers m and n. M=1 fictitious layers and N=18 order of expansion for the exponential matrix.
a
/
h
100
50
10
5
(
m
=
0
,
n
=
1
)
I mode
0.3318
0.3318
0.3312
0.3293
II mode
64.944
32.472
6.4944
3.2472
III mode
109.78
54.888
10.977
5.4870
(
m
=
1
,
n
=
0
)
I mode
2.9861
2.9846
2.9360
2.8012
II mode
194.83
97.417
19.483
9.7417
III mode
329.33
164.66
32.908
16.415
(
m
=
1
,
n
=
1
)
I mode
3.3179
3.3159
3.2562
3.0923
II mode
205.37
102.69
20.537
10.269
III mode
347.14
173.56
34.685
17.297
(
m
=
1,  n=2)
I mode
4.3130
4.3097
4.2099
3.9450
II mode
234.16
117.08
23.416
11.708
III mode
395.80
197.89
39.537
19.699
(
m
=
2
,
n
=
1
)
I mode
12.270
12.243
11.497
9.9094
II mode
395.04
197.52
39.504
19.752
III mode
667.72
333.83
66.559
32.908
(
m
=
2, n=2)
I mode
13.264
13.233
12.369
10.572
II mode
410.74
205.37
41.074
20.537
III mode
694.26
347.10
69.186
34.170Table 8
Benchmark 4. Simply supported orthotropic rectangular plate. First three vibration modes in terms of no-dimensional circular frequencyω-=ω(a2/h)ρ/E2 for several combinations of half-wave numbers m and n. M=1 fictitious layers and N=18 order of expansion for the exponential matrix.
a
/
h
100
50
10
5
100
50
10
5
(
m
=
0,  n=1),  θ=0°
(
m
=
0,  n=1),  θ=90°
I mode
0.3173
0.3173
0.3167
0.3148
1.1130
1.1126
1.0987
1.0586
II mode
75.922
37.961
7.5922
3.7961
75.922
37.961
7.5922
3.7961
III mode
104.98
52.483
10.495
5.2459
368.24
184.12
36.822
18.408
(
m
=
1,  n=0),  θ=0°
(
m
=
1,  n=0),  θ=90°
I mode
10.007
9.9712
9.0075
7.2059
2.8553
2.8537
2.8058
2.6731
II mode
227.77
113.88
22.777
11.388
227.77
113.68
22.777
11.388
III mode
1104.7
552.35
110.41
46.861
314.89
157.44
31.457
15.679
(
m
=
1,  n=1), θ=0°
(
m
=
1,  n=1),  θ=90°
I mode
10.127
10.091
9.1126
7.2879
3.4236
3.4213
3.3511
3.1611
II mode
249.70
124.85
24.967
12.480
312.77
156.38
31.251
15.586
III mode
1107.6
553.78
110.69
47.184
441.11
220.55
44.099
22.031
(
m
=
1,  n=2),  θ=0°
(
m
=
1,  n=2),  θ=90°
I mode
10.538
10.500
9.4774
7.5865
6.1047
6.0947
5.8032
5.1319
II mode
306.16
153.08
30.604
15.284
341.88
170.94
34.157
17.031
III mode
1116.1
558.04
111.54
48.130
774.35
387.17
77.407
38.652
(
m
=
2,  n=1),  θ=0°
(
m
=
2,  n=1),  θ=90°
I mode
40.001
39.445
28.902
18.694
11.868
11.841
11.093
9.5179
II mode
466.89
233.44
46.684
23.337
560.57
280.28
56.021
27.943
III mode
2210.8
1105.3
187.77
51.133
656.59
328.31
65.435
32.318
(
m
=
2,  n=2),  θ=0°
(
m
=
2,  n=2),  θ=90°
I mode
40.364
39.801
29.152
18.872
13.685
13.649
12.644
10.634
II mode
499.39
249.69
49.921
24.938
625.52
312.73
62.343
30.811
III mode
2215.1
1107.5
188.74
51.982
882.21
441.09
88.126
43.879Table 9
Benchmark 5. Simply supported isotropic cylindrical shell panel. First three vibration modes in terms of no-dimensional circular frequencyω-=ω(a2/h)ρ/E for several combinations of half-wave numbers m and n. M=100 fictitious layers and N=3 order of expansion for the exponential matrix.
R
α
/
h
1000
100
10
5
(
m
=
0
,
n
=
1
)
I mode
1067.1
106.71
10.683
5.3415
II mode
1068.3
106.83
10.690
5.3727
III mode
1855.7
185.57
18.557
9.2797
(
m
=
1
,
n
=
0
)
I mode
2.5478
2.5184
2.4853
2.3948
II mode
2040.3
204.03
20.411
10.216
III mode
3635.3
363.51
36.229
17.917
(
m
=
1
,
n
=
1
)
I mode
222.72
22.503
3.8758
3.2156
II mode
2364.5
236.45
23.656
11.842
III mode
4017.4
401.73
40.060
19.846
(
m
=
1
,
n
=
2
)
I mode
552.89
55.574
7.7851
5.7590
II mode
3031.5
303.16
30.328
15.181
III mode
5048.7
504.86
50.367
24.983
(
m
=
2
,
n
=
1
)
I mode
70.382
14.052
11.533
10.021
II mode
4230.8
423.08
42.325
21.174
III mode
7214.2
721.38
71.724
35.120
(
m
=
2
,
n
=
2
)
I mode
232.90
27.451
13.793
11.643
II mode
4639.2
463.92
46.414
23.226
III mode
7847.2
784.68
78.020
38.182Table 10
Benchmark 6. Simply supported orthotropic cylindrical shell panel. First three vibration modes in term of no-dimensional circular frequencyω-=ω(a2/h)ρ/E2 for several combinations of half-wave numbers m and n. M=100 fictitious layers and N=3 order of expansion for the exponential matrix.
R
α
/
h
1000
100
10
5
1000
100
10
5
(
m
=
0
,
n
=
1
)
,
θ
=
0
°
(
m
=
0
,
n
=
1
)
,
θ
=
9
0
°
I mode
1248.9
124.89
12.489
6.2444
1096.5
109.67
11.281
6.0108
II mode
1721.6
172.16
17.230
8.6360
1248.9
124.89
12.489
6.2444
III mode
3858.5
385.85
38.590
19.293
6057.8
605.78
60.570
30.272
(
m
=
1
,
n
=
0
)
,
θ
=
0
°
(
m
=
1
,
n
=
0
)
,
θ
=
9
0
°
I mode
8.4450
8.4392
7.6484
6.1872
2.3272
2.4079
2.3754
2.2863
II mode
2385.2
238.52
23.861
11.941
2385.2
238.52
23.861
11.943
III mode
12194
1219.4
121.48
51.358
3476.0
347.58
34.633
17.114
(
m
=
1
,
n
=
1
)
,
θ
=
0
°
(
m
=
1
,
n
=
1
)
,
θ
=
9
0
°
I mode
230.37
24.568
8.0953
6.3709
352.48
35.490
5.3604
4.1227
II mode
2939.9
294.00
29.454
14.782
3588.6
358.85
35.788
17.735
III mode
12257
1225.7
122.10
52.106
6559.0
655.90
65.574
32.756
(
m
=
1
,
n
=
2
)
,
θ
=
0
°
(
m
=
1
,
n
=
2
)
,
θ
=
9
0
°
I mode
641.70
64.887
10.877
6.8670
644.27
65.496
12.428
9.0851
II mode
4155.7
415.58
41.604
20.824
4164.3
416.43
41.609
20.709
III mode
12445
1244.4
123.94
54.240
12370
1237.0
123.62
54.652
(
m
=
2
,
n
=
1
)
,
θ
=
0
°
(
m
=
2
,
n
=
1
)
,
θ
=
9
0
°
I mode
79.744
39.112
28.352
18.587
160.90
20.150
11.545
9.9019
II mode
5062.8
506.29
50.681
25.371
6476.9
647.66
64.461
31.673
III mode
23492
2349.1
205.91
56.149
7983.5
798.33
79.684
39.571
(
m
=
2
,
n
=
2
)
,
θ
=
0
°
(
m
=
2
,
n
=
2
)
,
θ
=
9
0
°
I mode
243.62
46.084
29.034
19.077
364.43
40.611
16.527
13.080
II mode
5849.5
584.96
58.566
29.293
6951.1
695.08
69.132
33.884
III mode
23598
2359.6
208.24
58.169
13115
1311.5
131.01
62.003Table 11
Benchmark 7. Simply supported isotropic cylinder. First three vibration modes in term of no-dimensional circular frequencyω-=ω(a2/h)ρ/E for several combinations of half-wave numbers m and n. M=100 fictitious layers and N=3 order of expansion for the exponential matrix.
R
α
/
h
1000
100
10
5
(m=0, n=1)
I mode
7691.7
769.17
76.917
38.459
II mode
12342
1234.2
123.42
61.710
III mode
41586
4158.6
416.30
208.81
(
m
=
2, n=1)
I mode
2445.9
244.59
24.492
12.295
II mode
27270
2727.0
272.81
136.56
III mode
59189
5918.8
591.08
294.28
(
m
=
2,  n=2)
I mode
7712.1
771.23
77.252
38.820
II mode
33375
3337.5
333.85
167.09
III mode
61353
6135.3
612.74
305.15
(
m
=
2,  n=3)
I mode
13536
1353.7
135.69
68.323
II mode
40081
4008.1
400.89
200.56
III mode
65490
6548.9
654.17
325.96Table 12
Benchmark 8. Simply supported orthotropic cylinder. First three vibration modes in terms of no-dimensional circular frequencyω-=ω(a2/h)ρ/E2 for several combinations of half-wave numbers m and n. M=100 fictitious layers and N=3 order of expansion for the exponential matrix.
R
α
/
h
1000
100
10
5
1000
100
10
5
(
m
=
0,  n=1), θ=0°
(
m
=
0, n=1), θ=90°
I mode
8991.9
899.19
89.919
44.960
8991.9
899.19
89.919
44.959
II mode
12402
1240.2
124.16
62.281
39097
3909.7
391.39
196.33
III mode
138826
13883
1389.6
696.66
44039
4403.9
440.45
220.32
(
m
=
2, n=1), θ=0°
(
m
=
2, n=1), θ=90°
I mode
2486.4
248.65
24.921
12.541
5295.1
529.53
53.084
26.738
II mode
31697
3169.8
317.14
158.81
50390
5039.0
503.83
251.80
III mode
196438
19643
1961.3
974.72
58023
5802.3
579.62
288.86
(
m
=
2,  n=2),  θ=0°
(
m
=
2, n=2),  θ=90°
I mode
8071.9
807.21
80.903
40.715
11839
1184.0
119.10
60.562
II mode
38989
3898.9
390.16
195.44
56607
5660.7
565.47
281.81
III mode
196773
19677
1964.6
976.14
92402
9240.2
924.00
461.97
(
m
=
2, n=3),  θ=0°
(
m
=
2,  n=3), θ=90°
I mode
14640
1464.1
146.76
73.903
17656
1765.9
178.90
92.578
II mode
48229
4823.0
482.66
241.82
58769
5876.8
587.41
293.24
III mode
197336
19733
1970.0
978.55
134283
13428
1342.8
671.29Table 13
Benchmark 9. Simply supported isotropic spherical shell panel. First three vibration modes in terms of no-dimensional circular frequencyω-=ω(a2/h)ρ/E for several combinations of half-wave numbers m and n. M=100 fictitious layers and N=3 order of expansion for the exponential matrix.
R
α
/
h
1000
100
10
5
(
m
=
0, n=1)
I mode
999.01
99.924
10.213
5.3986
II mode
2040.3
204.02
20.360
10.115
III mode
3785.7
378.55
37.677
18.556
(
m
=
1,  n=0)
I mode
999.01
99.924
10.213
5.3986
II mode
2040.3
204.02
20.360
10.116
III mode
3785.7
378.55
37.677
18.556
(
m
=
1,  n=1)
I mode
1046.3
104.76
11.564
6.9554
II mode
2885.4
288.53
28.793
14.302
III mode
5111.6
511.13
50.821
24.937
(
m
=
1, n=2)
I mode
1076.3
108.52
16.906
12.400
II mode
4562.2
456.21
45.524
22.598
III mode
7857.9
785.72
77.933
37.847
(
m
=
2, n=1)
I mode
1076.3
108.52
16.906
12.400
II mode
4562.2
456.21
45.524
22.598
III mode
7857.9
785.72
77.933
37.847
(
m
=
2,  n=2)
I mode
1084.1
110.77
23.273
17.513
II mode
5770.8
577.07
57.581
28.565
III mode
9869.7
986.87
97.645
46.869Table 14
Benchmark 10. Simply supported orthotropic spherical shell panel. First three vibration modes in terms of no-dimensional circular frequencyω-=ω(a2/h)ρ/E2 for several combinations of half-wave numbers m and n. M=100 fictitious layers and N=3 order of expansion for the exponential matrix.
R
α
/
h
1000
100
10
5
1000
100
10
5
(
m
=
0, n=1), θ=0∘
(
m
=
0,  n=1), θ=90∘
I mode
2385.2
238.51
23.802
11.825
1037.9
104.08
12.585
7.7554
II mode
2902.0
290.20
29.026
14.513
2385.2
238.51
23.802
11.824
III mode
4371.6
437.15
43.559
21.518
12223
1222.3
121.65
51.888
(
m
=
1, n=0), θ=0∘
(
m
=
1, n=0), θ=90∘
I mode
1037.9
104.08
12.583
7.7554
2385.2
238.51
23.802
11.825
II mode
2385.2
238.51
23.802
11.824
2902.0
290.20
29.026
14.513
III mode
12223
1222.3
121.65
51.888
4371.6
437.15
43.559
21.518
(
m
=
1, n=1), θ=0∘
(
m
=
1, n=1), θ=90∘
I mode
1209.7
121.32
14.754
9.1755
1209.7
121.32
14.754
9.1755
II mode
4016.5
401.63
40.041
19.819
4016.5
401.63
40.041
19.819
III mode
12456
1245.5
123.88
54.506
12456
1245.5
123.88
54.506
(
m
=
1,  n=2), θ=0∘
(
m
=
1, n=2), θ=90∘
I mode
1682.9
169.01
21.986
14.305
1114.0
117.90
30.685
19.645
II mode
6838.7
683.84
68.035
33.373
5758.3
575.82
57.401
28.398
III mode
13147
1314.6
130.57
60.847
23599
2359.8
208.50
58.572
(
m
=
2,  n=1), θ=0∘
(
m
=
2, n=1), θ=90∘
I mode
1114.0
117.90
30.685
19.645
1682.9
169.01
21.986
14.305
II mode
5758.3
575.82
57.401
28.398
6838.7
683.84
68.035
33.373
III mode
23599
2359.8
208.50
58.572
13147
1314.6
130.57
60.847
(
m
=
2, n=2), θ=0∘
(
m
=
2, n=2), θ=90∘
I mode
1260.9
133.19
34.372
22.411
1260.9
133.19
34.372
22.411
II mode
8009.4
800.89
79.625
39.004
8009.4
800.89
79.625
39.004
III mode
23985
2398.3
216.62
65.349
23985
2398.3
216.62
65.349Figure 4
Benchmark 1, simply supported isotropic square plate. First three vibration modes in terms of displacement components through the thickness for thickness ratioa/h=10 and several combinations of half-wave numbers.Figure 5
Benchmark 2, simply supported orthotropic square plate(θ=0°). First three vibration modes in terms of displacement components through the thickness for thickness ratio a/h=10 and several combinations of half-wave numbers.Figure 6
Benchmark 2, simply supported orthotropic square plate (θ=90°). First three vibration modes in terms of displacement components through the thickness for thickness ratio a/h=10 and several combinations of half-wave numbers.Figure 7
Benchmark 7, simply supported isotropic cylinder. First three vibration modes in terms of displacement components through the thickness for thickness ratioRα/h=10 and several combinations of half-wave numbers.Figure 8
Benchmark 8, simply supported orthotropic cylinder(θ=0°). First three vibration modes in terms of displacement components through the thickness for thickness ratio Rα/h=10 and several combinations of half-wave numbers.Figure 9
Benchmark 8, simply supported orthotropic cylinder(θ=90°). First three vibration modes in terms of displacement components through the thickness for thickness ratio Rα/h=10 and several combinations of half-wave numbers.Figure 10
Benchmark 9, simply supported isotropic spherical shell panel. First three vibration modes in terms of displacement components through the thickness for thickness ratioRα/h=5 and several combinations of half-wave numbers.Figure 11
Benchmark 10, simply supported orthotropic spherical shell panel(θ=0°). First three vibration modes in terms of displacement components through the thickness for thickness ratio Rα/h=5 and several combinations of half-wave numbers.Figure 12
Benchmark 10, simply supported orthotropic spherical shell panel(θ=90°). First three vibration modes in terms of displacement components through the thickness for thickness ratio Rα/h=5 and several combinations of half-wave numbers.Table5 presents a square isotropic plate (benchmark 1). It shows the first three vibration modes for each possible combination of half-wave numbers (m,n) 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 (a/h=10) in terms of the three displacement components u*, v*, and w* 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 w* through the thickness is zero. The first mode has a constant transverse displacement w* through the thickness and linear in-plane displacements u* 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 n=0. The third vibration mode always has linear transverse displacement w* through the thickness for each pair of half-wave numbers. In-plane displacements u* 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 n.Table6 shows a square orthotropic plate (benchmark 2) and lists the first three vibration modes for each possible combination of half-wave numbers (m,n) 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 u*, v*, and w* for half-wave numbers (1,0), (1,1), (2,1), and (2,2) and for a thick plate (a/h=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 half-wave 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.Tables7 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 a and b are not the same (b=3a). 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).Table9 gives the first three vibration modes of the isotropic cylindrical shell panel (benchmark 5) for several combinations of half-wave numbers (m,n). The shell is isotropic but the frequencies are always different for each couple of (m,n) because of the geometry with different a and b dimensions. The same considerations are also valid for benchmark 6 in Table 10 (the orthotropic cylindrical shell panel). The geometry with different a and b 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.Tables11 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 (m,n) always has constant displacements u* and w*, and a quasi-constant v* displacement, but in the case m=0 and n=1 the only component different from zero is u*. The second mode always has constant displacements v* and w* and a linear u* displacement, but in the case m=0 and n=1u* is zero and w* and v* are constant. The third mode has quasi-constant u*, v*, and w* displacements, but in the case m=0 and n=1 the component u* 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 (m,n) 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 w* 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∘).Tables13 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 a and b. 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 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 theHα effect in cylindrical shell panels and both Hα and Hβ 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 Hα=1 (which means z/Rα≃0) for the benchmark 5 concerning the isotropic cylindrical shell panel. The approximation Hα=1 is valid for very thin shells (Rα/h equals 1000 or 100) even if it could have some problems for particular vibration modes and higher values of half-wave numbers (m,n). For thick shells the use of the exact geometry without the approximation Hα=1 is mandatory. Table 16 shows the comparison between the present 3D solution using the exact geometry and the same solution using the approximation Hα=Hβ=1 (which means z/Rα=z/Rβ≃0) for the benchmark 10 concerning the orthotropic spherical shell panel. The approximation Hα=Hβ=1 is valid for thin shells (Rα/h 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 z/Rα=z/Rβ≃0 is not valid.Table 15
Effects ofHα=(1+z/Rα) terms in benchmark 5 for simply supported isotropic cylindrical shell panel. First three vibration modes in terms of no-dimensional circular frequency ω-=ω(a2/h)ρ/E for several combinations of half-wave numbers m and n. M=100 fictitious layers and N=3 order of expansion for the exponential matrix.
R
α
/
h
1000
100
10
5
(
m
=
0, n=1)
I mode
1067.1
106.71
10.683
5.3415
II mode
1068.3
106.83
10.690
5.3727
III mode
1855.7
185.57
18.557
9.2797
I mode (Hα=1)
1067.1
106.71
10.686
5.3486
II mode (Hα=1)
1068.3
106.83
10.688
5.3684
III mode (Hα=1)
1855.7
185.56
18.550
9.2651
(
m
=
1,  n=0)
I mode
2.5478
2.5184
2.4853
2.3948
II mode
2040.3
204.03
20.411
10.216
III mode
3635.3
363.51
36.229
17.917
I mode (Hα=1)
2.7747
2.7473
2.7029
2.5828
II mode (Hα=1)
2040.3
204.03
20.403
10.201
III mode (Hα=1)
3635.3
363.52
36.303
18.076
(
m
=
1, n=1)
I mode
222.72
22.503
3.8758
3.2156
II mode
2364.5
236.45
23.656
11.842
III mode
4017.4
401.73
40.060
19.846
I mode (Hα=1)
222.73
22.535
4.0440
3.3738
II mode (Hα=1)
2364.5
236.45
23.668
11.868
III mode (Hα=1)
4017.4
401.73
40.115
19.965Table 16
Effects ofHα=(1+z/Rα) and Hβ=(1+z/Rβ) terms in benchmark 10 for simply supported orthotropic spherical shell panel. First three vibration modes in term of no-dimensional circular frequency ω-=ω(a2/h)ρ/E2 for several combinations of half-wave numbers m and n. M=100 fictitious layers and N=3 order of expansion for the exponential matrix.
R
α
/
h
1000
100
10
5
1000
100
10
5
(
m
=
0, n=1), θ=0°
(
m
=
0, n=1), θ=90°
I mode
2385.2
238.51
23.802
11.825
1037.9
104.08
12.585
7.7554
II mode
2902.0
290.20
29.026
14.513
2385.2
238.51
23.802
11.824
III mode
4371.6
437.15
43.559
21.518
12223
1222.3
121.65
51.888
I mode (Hα=Hβ=1)
2385.2
238.52
23.858
11.939
1037.9
104.09
12.602
7.7410
II mode (Hα=Hβ=1)
2902.0
290.20
29.030
14.529
2385.2
238.52
23.856
11.934
III mode (Hα=Hβ=1)
4371.6
437.15
43.623
21.667
12223
1222.3
122.04
51.743
(
m
=
1, n=0), θ=0°
(
m
=
1, n=0), θ=90°
I mode
1037.9
104.08
12.583
7.7554
2385.2
238.51
23.802
11.825
II mode
2385.2
238.51
23.802
11.824
2902.0
290.20
29.026
14.513
III mode
12223
1222.3
121.65
51.888
4371.6
437.15
43.559
21.518
I mode (Hα=Hβ=1)
1037.9
104.09
12.602
7.7410
2385.2
238.52
23.858
11.939
II mode (Hα=Hβ=1)
2385.2
238.52
23.856
11.934
2902.0
290.20
29.030
14.529
III mode (Hα=Hβ=1)
12223
1222.3
122.04
51.743
4371.6
437.15
43.623
21.667
(
m
=
1, n=1), θ=0°
(
m
=
1, n=1), θ=90°
I mode
1209.7
121.32
14.754
9.1755
1209.7
121.32
14.754
9.1755
II mode
4016.5
401.63
40.041
19.819
4016.5
401.63
40.041
19.819
III mode
12456
1245.5
123.88
54.506
12456
1245.5
123.88
54.506
I mode (Hα=Hβ=1)
1209.7
121.33
14.810
9.2169
1209.7
121.33
14.810
9.2169
II mode (Hα=Hβ=1)
4016.5
401.64
40.137
20.022
4016.5
401.64
40.137
20.022
III mode (Hα=Hβ=1)
12456
1245.6
124.29
54.312
12456
1245.6
124.29
54.312


## 5. 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. 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.



---

*Source: 479738-2014-03-09.xml*