Free Vibration Analysis for Shells of Revolution Using an Exact Dynamic Stiffness Method

An exact generalised formulation for the free vibration of shells of revolutionwith general shapedmeridians and arbitrary boundary conditions is introduced. Starting from the basic shell theories, the vibration governing equations are obtained in the Hamilton form, from which dynamic stiffness is computed using the ordinary differential equations solver COLSYS. Natural frequencies and modes are determined by employing the Wittrick-Williams (W-W) algorithm in conjunction with the recursive Newton’s method, thus expanding the applications of the abovementioned techniques from one-dimensional skeletal structures to two-dimensional shells of revolution. A solution for solving the number of clamped-end frequencies J 0 in the W-W algorithm is presented for both uniform and nonuniform shell segment members. Based on these theories, a FORTRAN program is written. Numerical examples on circular cylindrical shells, hyperboloidal cooling tower shells, and spherical shells are given, and error analysis is performed. The convergence of the proposed method on J 0 is verified, and comparisons with frequencies from existing literature show that the dynamic stiffness method is robust, reliable, and accurate.


Introduction
Shells of revolution form an important class of structures that are used in a variety of engineering applications, for example, pipes, chimneys, cooling towers, containment vessels, and aircraft fuselages.These structures usually operate in complex conditions subject to dynamic loads.In order to understand the dynamic performance of shells of revolution better, exact modal analysis is necessary.
From a historical perspective, shell theories can be traced back to the 19th century when Clapeyron and Lamé [1] formulated the membrane and bending theories of shells.Since Love [2] established the basis of thin elastic shell theory, much effort was expended on its development and attempts have been made to investigate the free vibration of shells.Arnold and Warburton [3] solved the flexural vibration of thin cylindrical shells with simple boundary conditions based on Rayleigh's principle and showed that the lowest natural frequency is usually obtained with a relatively high circumferential wave number.Forsberg [4] studied comprehensively the influences of boundary conditions on the modal characteristics of thin cylindrical shells by using the Fourier series for the displacement functions based on Flügge's theory.
Leissa [5] systematically reviewed and summarised the vibration of different shell theories and shell types in his monograph.Luah and Fan [6] analysed the free vibration of shells of revolution using the spline finite element method.Tan [7] introduced a substructuring method for predicting the natural frequencies of general shells of revolution.Gautham and Ganesan [8] investigated the free vibration characteristics of isotropic spherical caps by incorporating a first-order shear deformable semianalytical shell finite element.Bardell et al. [9,10] applied an h-p version of the finite element method to the free vibration of open cylindrical and conical isotropic shell panels.Mantari et al. [11] performed free vibration analysis of multilayered shells by using a new accurate higher-order shear deformation theory, and a three-dimensional solution was presented by Jin et al. [12] on free vibration of isotropic conical shells with elastic boundary restraints.Although the free vibration of isotropic shells of revolution has been extensively examined, there is still continuing relevant research conducted with new methods [13,14].
Besides simple shaped and isotropic shells, free vibration of doubly curved shells becomes popular and has been investigated extensively in recent years.Tornabene [15] studied the dynamic behaviour of laminated composite doubly curved shells of revolution using a Generalised Differential Quadrature (GDQ) method, and very good agreement was observed by comparing with results from literature and commercial FEM software.Later, the GDQ method was applied to the free vibration of laminated and anisotropic doubly curved shells of revolution with a free-form meridian [16][17][18].In Bacciocchi et al. [19] and Tornabene et al. [20], the GDQ method was further extended to the computation of natural frequencies of doubly curved shells with variable thickness.
This paper employs an exact dynamic stiffness method and the Wittrick-Williams (W-W) algorithm [21,22].It considers structures as continuous systems with infinite degrees of freedom and hence dynamic stiffness is obtained from exact vibration governing equations.Natural frequencies and modes are solved by using the W-W algorithm in conjunction with the recursive Newton's method [23,24].Since governing equations are solved exactly when calculating the dynamic stiffness, the proposed method is exact.
The dynamic stiffness method and the W-W algorithm are widely used in free vibration analysis of skeletal structures, for example, nonuniform Timoshenko beams [28], three-layered sandwich beams [29], and functionally graded beams [30].However, there are few applications on plates and shells.A recent paper [31] applied the dynamic stiffness method and the W-W algorithm to the free vibration analysis of thin circular cylindrical shells.This work, which is based on a thesis of the first author [32], extends the research of [31] substantially.A generalised formulation for free vibration of shells of revolution with general shaped meridians and arbitrary boundary conditions is presented, enriching the literature on shell vibration using the dynamic stiffness method.Numerical results on the vibration behaviour of circular cylindrical shells, cooling tower shells, and spherical shells are given, showing that the exact dynamic stiffness method is applicable, accurate, and robust.

Basic Equations
A thin, homogeneous, isotropic, and circumferentially closed shell of revolution is considered in this paper.Based on the Kirchhoff-Love hypothesis, the vibration of the shell is represented by the middle surface displacements (u, v, w) along the meridian, circumferential, and radial directions, respectively.
Referring to Figure 1,  is the angle between the normal direction ⃗  of the meridian at an arbitrary point  and the positive direction of -axis;  is the angle between the meridian plane at point  and a base meridian plane OXZ;  1 and  2 are the centres of the meridian curvature and the parallel circle; r and  are the radii of the parallel circle and meridian curvature at point P, respectively.Young's modulus of the material is E, Poisson's ratio is ], density is , and the thickness is h.
According to the thin shell theory, the strain-displacement relationship of a shell of revolution is where the strain vector {} = {  ,   ,   ,   ,   ,   }  and the displacement vector {Δ} = {, V, }  .The expression of [L] can be found in [5].Similarly, the force-strain relationship of a shell of revolution is where the internal force vector {N} = {  ,   ,   ,   ,   ,   }  and [D] is the stiffness matrix given by Denote by dA the area of an infinitesimal portion in the middle surface of a shell of revolution (see Figure 2).The area of this infinitesimal portion is  = .The total strain energy  for free vibration of the shell of revolution can be integrated along the whole middle surface S: and the total free vibration kinetic energy is obtained as Assume that the dynamic displacement functions u, v, and  are where  is the circumferential wave number and  is a circular frequency.According to Hamilton's principle, Since time interval ( 1 ,  2 ) is arbitrary, the symplectic components (generalised displacements and forces) can be obtained from (7) after displacement substitution and integration manipulation.
Referring to Figure 3(a), when integrating along meridian ⃗  of a shell of revolution with a positive Gaussian curvature,  1 <  2 ,  > 0, and  > 0. The generalised displacements and forces obtained from (7) are given by where ()  denotes differentiation with respect to  and The governing equations of motion in free vibration in the Hamilton form are given by where , I is a fourth-order identity matrix, and [S] is a symmetric eighth-order matrix.For positive Gaussian curvature shells, [S] = [S P ] as is given in the Appendix.Similarly, for a shell of revolution with a negative Gaussian curvature (see Figure 3(b)),  1 >  2 ,  < 0, and  > 0. The only change in the symplectic components is a minus before the generalised forces.Equation ( 10) is still applicable should the matrix [S] be replaced by [S N ], which is also given in the Appendix.

Dynamic Stiffness Matrices
Should the circumferential wave number  take a specific value, ( 10) is one-dimensional with respect to  only.Therefore, a thin shell of revolution can be analysed as a general  one-dimensional skeleton-like structure with four degrees of freedom ( 1 ,  2 ,  3 ,  4 ) and four corresponding internal forces (V 1 , V 2 , V 3 , V 4 ).The dynamic stiffness matrix of free vibration of a shell of revolution can now be set up in a similar way as that in the classical one-dimensional skeletal theory.
Referring to Figure 4, a shell of revolution can thus be divided into several shell segments along its meridian.Taking the example of segment element (e), the angle coordinates at its two ends are   and   , and the segment-end displacement vector {d}  and force vector {F}  are defined as Mathematically apply the boundary conditions in ( 12) to (10) in turns: where {  } is a unit vector with the jth element equal to one, and the corresponding element dynamic stiffness matrix [k i ]  is formulated by Due to the complexity, ( 10) is solved numerically using adaptive ordinary differential equations (ODEs) solver COLSYS [33,34].The global dynamic stiffness matrix K of a shell of revolution is assembled in a direct way with regular element location vectors.

The Wittrick-Williams Algorithm
The W-W algorithm gives the number of frequencies below a specific value of trial frequency  * of interest with the expression of where  is the total number of frequencies exceeded by a specific trial frequency  * ,   = {K( * )} is the sign count of the global dynamic stiffness matrix K and equals the number of negative elements on the diagonal of an upper triangular matrix obtained from K by applying standard Gaussian elimination without row exchange, and  0 is the number of clamped-end frequencies exceeded by  * and can be accumulated from   0 of each shell segment element over the whole meridian as is expressed in In the present work, a substructuring method is employed by taking advantage of the self-adaptability of COLSYS.When computing the dynamic stiffness matrix of a shell segment element, for example, element (e), couple of submeshes (termed as (ê)) are defined by COLSYS as is shown in Figure 5.
Applying the W-W algorithm again on the submesh (ê), Since COLSYS is capable of controlling the error tolerance due to its self-adaptability, the number of clamped-end frequencies exceeded by  * on submesh (ê) has to be zero, suggesting that  (ê) 0 = 0. Otherwise, the clamped-end modes can be enlarged infinitely, resulting in the unsatisfaction of the given error tolerance.
Take a submesh element (α  , α+1 ) on element (e), for instance; its corresponding dynamic stiffness matrix [ ki ] (ê) can be obtained by linearly combining the eight solutions of ( 10) and (12).Note that the end displacement and force vectors for submesh are {d} (ê) and {F} (ê) , respectively; the submesh element stiffness matrix [ ki ] (ê) satisfies By applying the boundary conditions in (12) in turns again, ( 17) is expanded as where As  (ê) 0 = 0, [B] (ê) is nonsingular (otherwise  * is one of the clamped-end frequencies, which is conflicted with  (ê) 0 = 0), and ( 20) is obtained: Since K (ê) ( * ) can be obtained from regular assembly of [ ki ] (ê) , the number of clamped-end frequencies of the shell segment (e) can now be computed from (16).For a general shell of revolution, each shell segment element is different and the proposed method is capable of solving  0 for nonuniform members.It can be verified that, with the decrease of the length of a shell segment along the meridian,  0 converges to zero.
By employing the W-W algorithm repeatedly, the upper and lower bound of each frequency under a given circumferential wave number  can be predicted.Recursive Newton's method is employed to obtain the exact natural frequencies and modes under the specific n, with its details in [23,24].Based on these theories, a FORTRAN program was written by the authors to determine the eigenproblems of a shell of revolution with general shaped meridians and arbitrary boundary conditions, and numerical examples provided in this paper are all computed by this program.

Numerical Examples
This section presents applications of free vibration analysis on circular cylindrical shells, hyperboloidal shells, and spherical shells using the proposed dynamic stiffness method. 0 ,   , and  are computed, showing the convergence of the method presented in the previous section.

Circular Cylindrical Shells.
As the formulae are derived from general shells of revolution, the proposed dynamic stiffness method can be simply employed for the free vibration analysis of circular cylindrical shells.The infinitesimal distance along the meridian is  =  with  = 90 ∘ and  = ∞.The radius  of parallel circles is constant and equals the radius of the circular cylindrical shell, and the governing equations can be written in the Hamilton form of  divided equally along the meridian. 0 ,   , and  with different number of segment members ne at  = 2,  * = 3000 rad/s are given in Table 1.It can be observed that the sum of  0 and   always converges to .With the increase of ne,  0 will be reduced to zero, suggesting that the calculation of  0 can be avoided should each shell segment member be short enough.In the following studies, ne is always large enough so that  0 = 0 can be guaranteed if not specifically stated.

Axisymmetric Vibrations (𝑛 = 0). For axisymmetric vibration of a Shear Diaphragm-Shear Diaphragm (SD-SD) circular cylindrical shell whose boundary conditions are
where   and   are the moment and force along the meridian, Leissa [5] gave the analytical solutions for both the pure-torsion and the axial-radial coupled vibrations.Table 2 lists the lowest twenty natural frequencies obtained from both the present method and Leissa [5], and very good agreement is reached.It is worth mentioning that the dynamic stiffness approach identified the zero frequency.To investigate the errors between the present method and the results from existing literature, a notation  in ( 23) is introduced: where Ω current and Ω literature are the frequencies obtained from the current method and the existing literature, respectively.Consequently, log 10 () is the error index.Errors for those frequencies which are exactly the same as the values in the literature due to limited significant digits are not computed in this way, but adopting (−d) instead, where  is the number of significant digits for the presented frequencies.Referring to Figure 6, which shows the discrete error points between the current dynamic stiffness method and Leissa [5] without the consideration of zero frequency, conclusions arrive such that the precision of the proposed method is very high, with errors smaller than the level of 10 −3 comparing with the analytical solutions.It is also observed that the dynamic stiffness method produces better results for axialradial coupled frequencies, as all the points in Figure 6 with  > 10 −6 belong to the pure-torsion vibrations.Denote by  the number of half-waves in the axial direction.It is known that, for each m, there are two vibration modes: pure-torsion and axial-radial coupled vibration.Figure 7 shows the first three pure-torsion and axialradial coupled modes of the axisymmetric vibration with the exclusion of zero frequency.The modes are normalised to the maximum displacement and the trends are as expected.It is apparent from Figure 7 that although the amplitude of radial displacement  increases with the increase of m, displacement  along the axial direction is still dominant for those axial-radial coupled vibrations when  is small.Chung [25] analysed the free vibration of a circular cylindrical shell with clamped-free (C-F) boundary condition using a general analytical method based on Sander's shell theory.
The geometrical parameters of the circular cylindrical shell are as follows: length L = 0.5112 m, radius  = 0.2162 m, and thickness h = 0.0015 m; the material properties are as follows: Young's modulus  = 1.83 × 10 11 N/m 2 , Poisson's ratio ] = 0.3, and density  = 7492 kg/m 3 .The first three natural frequencies with  from 0 to 4 are tabulated and compared with results from Chung [25] in Table 3, reaching very good agreement for most data points.Figure 8 gives the first three normalised vibration modes at  = 2, which are radial-dominant and the maximum radial displacements are always reached at the free end.Together with the SD-SD shell examined previously, good agreement and small errors are achieved in both cases, suggesting that the dynamic stiffness method is insensitive to the sizes or boundary conditions of shells.

Cooling Tower Shells.
Hyperbolic meridians are commonly used in cooling tower shells, which are typical structures with negative Gaussian curvatures.The base of a cooling tower is usually well attached to the ground so that the boundary condition is clamped-free.The geometry of a cooling tower is schematically shown in Figure 9 with its meridian equation given as   Table 4 presents  0 ,   , and  of the cooling tower shell with different ne at  = 0,  * = 100 rad/s, showing that the present method is well capable of capturing  0 for nonuniform shell segments.
Table 5 lists the first three natural frequencies of the cooling tower shell with  = 0, 1, and 2. The tower was also analysed by Deb Nath [26] using a ring finite element approach based on Vlasov and Love's constitutive laws, with   results given opposite in Table 5.It is observed that results from the present method agree well with that from Deb Nath [26], demonstrating that the dynamic stiffness method is capable of dealing with vibrations of negative Gaussian curvature shells of revolution.

Spherical Shells.
For spherical shells, R is constant and held positive, and  = sin .de Souza and Croll [27] investigated two apex-closed spherical shell caps with their bottom clamped (see Figure 10) using an energy variation procedure.To avoid singularity of [L] at the apex ( = 0), a small circular hole   = 0.125%   is introduced by de Souza and Croll [27].The same holes are created in the present dynamic stiffness analysis.The first four nondimensional natural frequencies (in  0 = √/ 2 (1 − ] 2 )) with  varying from 1 to 3 of the two spherical shell caps are given in Tables 6 and 7, respectively.
The data from the present method are in satisfactory agreement with those from de Souza and Croll [27], showing that the present method is applicable and accurate for the free vibration analysis of positive Gaussian curvature shells, for example, spherical shells.
Error analysis on the spherical shells in Figure 11 shows that generally a 10 −2 magnitude error is achieved by the proposed dynamic stiffness method.Further investigations on Figure 11 demonstrate that (1) the proposed dynamic stiffness approach is not sensitive to the shape of shells as similar levels of errors are achieved for both the 45 ∘ and 90 ∘ spherical shell caps and (2) errors do not necessarily increase along with the circumferential wave number  or half-wave number m, suggesting that the accuracy of higherorder frequencies is not affected by the presented dynamic stiffness method.

Conclusions
An exact dynamic stiffness approach is formulated for the free vibration of shells of revolution in this reported work.
The vibration equations are degraded into a series of onedimensional problems with respect to different circumferential wave number n, thus expanding the applications of the W-W algorithm from one-dimensional skeletal structures to two-dimensional shells of revolution.Since the governing equations are solved exactly by ODEs solver COLSYS when computing the dynamic stiffness matrices, the proposed method is exact.A solution for calculating the number of clamped-end frequencies  0 in the W-W algorithm is also proposed and its convergence is validated.
Based on these theories, a FORTRAN program is written.Numerical examples on free vibration of circular cylindrical shells, cooling tower shells, and spherical shells cover the structures with both positive and negative Gaussian curvatures, as well as straight meridians.Conclusions are reached as follows: (a) The proposed dynamic stiffness method which employs a generalised formulation is applicable to a variety of shells of revolution.(b) Error analysis on results computed from the dynamic stiffness method shows that the present method is satisfactory (generally < 10 −2 ), providing an accurate way of estimating the engineering eigenproblem of shells of revolution.(c) The dynamic stiffness method is insensitive to the shape or the size of a shell, and the precision does not decrease for higher-order frequencies.Zero frequencies are also identified for SD-SD circular cylindrical shells.These merits show that the proposed method is robust for free vibration analysis.
In summary, the reported work in this paper demonstrates that the dynamic stiffness method is applicable, accurate, and robust for the free vibration analysis of shells of revolution.Though based on isotropic thin shell theories, the proposed method can be straight-forwardly extended to the free vibration of composite shells or shells with the consideration of shear deformation and rotary inertia, providing a wider perspective to future research and application in engineering.

Appendix
Consider the following:

Figure 1 :
Figure 1: Geometry of a shell of revolution.

2 Figure 2 :
Figure 2: An infinitesimal portion of a shell of revolution.

Figure 3 :
Figure 3: Change of the angle  according to different Gaussian curvatures.

Figure 6 :
Figure6: Error index of an axisymmetric vibration circular cylindrical shell between the present method and Leissa[5].

Figure 7 :
Figure 7: The first three pure-torsion and axial-radial coupled vibration modes at  = 0.

Figure 8 :
Figure 8: Normalised modes of the lowest three frequencies of the C-F circular cylindrical shell at  = 2.

Figure 11 :
Figure 11: Error index of spherical shell caps between the present method and de Souza and Croll [27].

Table 2 :
The lowest twenty natural frequencies (rad/s) of axisymmetric vibration.

Table 3 :
The first three natural frequencies (Hz) of a C-F circular cylindrical shell with n from 0 to 4.

Table 5 :
Natural frequencies (Hz) of a cooling tower shell.