# An Exact Dynamic Stiffness Formulation for Predicting Natural Frequencies of Moderately Thick Shells of Revolution

Authors: Xudong Chen; Kangsheng Ye
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

An exact dynamic stiffness formulation is proposed for calculating the natural frequencies of shells of revolution based on the first-order Reissner-Mindlin theory. Equations of motion are reduced to be one-dimensional, should the circumferential wave number is specified, and then are rewritten in Hamilton form. Therefore, a shell of revolution with a moderate thickness is analysed as a special skeletal structure with five degrees of freedom at element ends. Dynamic stiffnesses are constructed directly from the one-dimensional governing vibration equations using classic skeletal theory. Number of natural frequencies below a specific value is counted by applying the Wittrick-Williams (W-W) algorithm, and exact eigenvalues can be simply obtained using bisection method. A solution to the number of clamped-end frequenciesJ0 in the W-W algorithm is also proposed and proven to be reliable. Numerical examples on a variety of shells with different boundary conditions are investigated, and results are well compared and validated. Influences of a variety of shell parameters on natural frequencies of both cylindrical and spherical shells are discussed in detail, demonstrating that the proposed dynamic stiffness formulation is applicable to analyse natural frequencies of moderately thick shells of revolution with high accuracy.


---

## Body

## 1. Introduction

Among all the structure types, shells of revolution have been widely used in tubes, containment vessels, aircraft fuselages, etc. Their natural frequencies are important parameters when performing further dynamic structural analysis. Therefore, exact free vibration analysis of shells of revolution is of considerable scientific and engineering significance.As the simplest shell of revolution, circular cylindrical shells are the most widely employed in engineering applications as they are easy to manufacture. Since Love [1] established the basis of thin elastic shell theory, loads of effort have been devoted to the free vibration analysis of circular cylindrical shells over decades [2, 3]. Besides circular cylindrical shells, other types of shells like conical, spherical, elliptical, parabolic, and hyperbolic shells are also frequently used in civil, mechanical, aeronautical, and marine engineering. The earliest study on free vibration of elliptical shells can be traced back to Lamb [4] in 19th century. Following researchers [5, 6] investigated the vibration behaviour of both prolate and oblate spheroidal shells. Lin and Lee [7] and Hoppmann et al. [8] studied the free vibration of paraboloidal shells. Carter et al. [9] presented a finite element approach to calculate the natural frequencies of hyperbolic shells of revolution. Li et al. [10] solved the natural frequencies of conical shells with Hamilton’s principle and Rayleigh-Ritz method. In recent decades, general solutions to free vibration of shells of revolution were provided; e.g., Luah and Fan [11] analysed the free vibration of shells of revolution by introducing the spline finite element method; Tan [12] predicted the natural frequencies of general shells of revolution using a substructuring approach.However, all the studies above are based on thin shell theories. For shells with moderate thickness, the employment of Reissner-Mindlin assumptions with the consideration of transverse shear deformation and rotary inertia is necessary. There are some researches conducted with moderately thick shell theories. Naghdi and Cooper [13] analysed the free vibration of circular cylindrical shells including shear deformation and rotary inertia using an elastic wave propagation approach. Mirsky and Herrmann [14] investigated nonaxially symmetric vibration motions of moderately thick circular cylindrical shells. Relevant research was collected by Leissa in his monograph [15].Recently, Hosseini-Hashemi and Fadaee [16] obtained the closed-form solutions to the free vibration of moderately thick spherical shell panels using new auxiliary and potential functions in conjunction with the separation method of variables. Su et al. [17] analysed the free vibration of moderately thick functionally graded open shells with general boundary conditions using the Rayleigh-Ritz method and Fourier series. Tornabene and Viola [18] conducted free vibration examination on parabolic shells with Generalized Differential Quadrature (GDQ) method based on the first-order shear deformation theory. For even thicker shells, Kang and Leissa [19] employed three-dimensional theory to analyse the free vibration of thick hyperbolic shells with Ritz method. Al-Khatib and Buchanan [20] investigated the free vibration of a parabolic shell of revolution using the three-dimensional finite element method.Unlike thin shells, free vibration of moderately thick shells of revolution with general meridian shapes and arbitrary boundary conditions is not that extensively investigated due to its complexity in governing equations. In this paper, a dynamic stiffness approach for the free vibration of moderately thick shells of revolution is formulated in conjunction with the Wittrick-Williams (W-W) algorithm [21, 22]. It treats structures as continuous systems with infinite degrees of freedom and dynamic stiffnesses are calculated directly from vibration governing equations. Natural frequencies are solved computationally by employing the W-W algorithm together with the bisection method. Since the governing partial differential equations (PDE) are degraded analytically into ordinary differential equations (ODE) when calculating the dynamic stiffnesses, the proposed formulation is more accurate over traditional finite element approaches where the differential equations are discretized into linear algebraic equations.The dynamic stiffness method and the W-W algorithm have been used in free vibration analysis of skeletal structures such as nonuniform Timoshenko beams [23], three-layered sandwich beams [24] and functionally graded beams [25]. However, limited attention is paid to the plates and shells from the computational community. Pagani et al. [26] analysed the free vibration of thin-walled structures with exact dynamic stiffness elements. El-Kaabazi and Kennedy [27] applied the dynamic stiffness method and the W-W algorithm to calculate the frequencies and modes of thin circular cylindrical shells with variable thickness.The present authors proposed the idea of dynamic stiffness method in investigating the free vibration of shells of revolution [28], which was summarised in [29] based on the classic thin shell theory. Meanwhile, endeavours on predicting the free vibration characteristics of moderately thick circular cylindrical shells [30] and elliptical shells [31] were made. A further comparison study [32] on the dynamic stiffness approach between thin and moderately thick circular cylindrical shells was conducted, showing that frequencies from moderately thick circular cylindrical shells can be used as alternatives to those from thin shells with no observation of shear locking when the thickness to radius ratio is small.In this paper, a concise but complete formulation for free vibration of moderately thick shells of revolution with general shaped meridians and arbitrary boundary conditions is presented. Since the number of circumferential half-waves is an integer, the original vibration partial differential equations are degraded into series of ordinary differential equations, and hence W-W algorithm is employed. A solution ofJ0 in the W-W algorithm is provided. Numerical results on natural frequencies of moderately thick circular cylindrical shells, elliptical shells, and parabolic and hyperbolic shells are given, showing that the proposed dynamic stiffness formulation is applicable and accurate.


## 2. Equations of Motion

Consider a homogeneous, isotropic, and circumferentially closed shell of revolution with moderate thickness (see Figure1). Referring to Figure 1,u,v, andw are the translational displacements along the meridian, circumferential, and radial directions on the middle surface at pointP.φ is the angle between the normal at pointP andz-axis, andθ is the angle along the circumferential direction.r andR are the radii of the parallel circle and meridian curvature, respectively. Young's modulus isE, Poisson's ratio isν, density isρ, and the thickness ish.Figure 1
A general shell of revolution and the coordination.Based on the Reissner-Mindlin hypothesis, vibration of a moderately thick shell is represented by the middle surface displacement vectoru,v,w,ψφ,ψθT, where ψφ and ψθ are the angular displacements of the middle surface alongφ andθ directions, respectively.The strain-displacement relationship of a shell of revolution is given in(1)ε=LΔwhere the strain vector ε=εφ,εθ,εφθ,χφ,χθ,χφθ,γφ,γθT, and the displacement vector Δ=u,v,w,ψφ,ψθT. The differential matrix [L] can be referred to in [15].Similarly, the force-strain relationship is(2)N=Dεwhere the internal force vector N=Nφ,Nθ,Nφθ,Mφ,Mθ,Mφθ,Qφ,QθT and [D] is the stiffness matrix(3)D=D100D2where(4)D1=Eh1-ν21ν0000100001-ν2000h21200sym.h2120h21-ν24(5)D2=κ-Gh00κ-Gh.In (5), κ-=π2/12,G is the shear modulus, and G=E/2(1+ν).Referring to Figure2,dA denotes the area of an infinitesimal portion in the middle surface of a moderately thick shell of revolution, anddA=rRdφdθ. The total strain energyU for free vibration of a moderately thick shell can be obtained by integrating the middle surfaceS as(6)U=∬S12εTDεdAand the corresponding total kinetic energy is(7)T=ρ2∬Ωh∂u∂t2+∂v∂t2+∂w∂t2+h312∂ψx∂t2+∂ψθ∂t2dA.Hamilton’s principle suggests the total strain energy and kinetic energy satisfy(8)δ∫t1t2U-Tdt=0where time interval (t1,t2) is arbitrary. Displacement functionsu,v,w, ψφ, and ψθ can be assumed as(9)u=unφcos⁡nθsin⁡ωtv=vnφsin⁡nθsin⁡ωtw=wnφcos⁡nθsin⁡ωtψφ=ψφnφcos⁡nθsin⁡ωtψθ=ψθnφsin⁡nθsin⁡ωtn=0,1,2…wheren is the circumferential wave number andω is a circular frequency. Equation (9) is widely accepted and suitable for general vibration analysis of moderately thick shells of revolution. The Symplectic components (generalized displacements and forces) can be obtained through displacement substitution and integration manipulation. Similar to [29], when integrating along the meridian of a shell with a positive Gaussian curvature, generalized displacements and forces obtained from (8) are given by(10)u=u1=unu2=vnu3=wnu4=ψφnu5=ψθn(11)v=v1=Ehπ1-ν2νcos⁡φun+rRun′+νnvn+νsin⁡φ+rRwnv2=Ehπ1+ν-n2+h224nsin⁡φrRun-cos⁡φ2+h224sin⁡φcos⁡φrRvn+r2R+h224sin2⁡φrRvn′-h224nsin⁡φrψφn-h224sin⁡φcos⁡φrψθn+h224sin⁡φRψθn′v3=κ-Ghπ-rRun+rRwn′+rψφnv4=Ehπ1-ν2h212νcos⁡φψφn+rRψφn′+νnψθnv5=Ehπ1+νh224-nRun-cos⁡φRvn+sin⁡φRvn′-nψφn-cos⁡φψθn+rRψθn′where ′ denotes the differentiation with respect toφ. The free vibration governing equations can be obtained in the Hamilton form as(12)Jz′=Szwhere z=uTvTT; J=0-II0,I is a fifth-order identity matrix; [S] is a symmetric tenth-order matrix. For shells of revolution with a positive Gaussian curvature, [S] = SP. For those with a negative Gaussian curvature, [S] = [SN]. Details on SP and [SN] are given in the Appendix.Figure 2
An infinitesimal portion in a moderately thick shell of revolution.


## 3. Dynamic Stiffness

Governing vibration equations in (12) is one-dimensional with respect toφ when circumferential wave numbern is known. Consequently, generalized one-dimensional skeletal theory can be seamlessly extended to the two-dimensional moderately thick shells of revolution. Therefore, a moderately thick shell of revolution is treated as a special ‘beam’ with five degrees of freedom u1u2u3u4u5T and five corresponding internal force components v1v2v3v4v5T.Figure3 illustrates a possible mesh division of a moderately thick shell of revolution, which is divided into couple of shell elements along its meridian. Taking the instance of element (e), the element-end displacement vector de and force vector Fe are(13)de=u1φa,u2φa,u3φa,u4φa,u5φa,u1φb,u2φb,u3φb,u4φb,u5φbTand(14)Fe=-v1φa,-v2φa,-v3φa,-v4φa,-v5φa,v1φb,v2φb,v3φb,v4φb,v5φbTwhere φa and φb are angle coordinates at the element ends shown in Figure 3, respectively.Figure 3
Mesh divisions along the meridian.Should boundary conditions in (15) be applied to (12) in turn(15)de=ej,j=1,2…,10(ej is a unit vector with thej-th element equals to one), the corresponding element dynamic stiffness matrix [ki]e is obtained as(16)kie=Fj=1eFj=2eFj=3eFj=4eFj=5eFj=6eFj=7eFj=8eFj=9eFj=10e.Equation (12) is too complicated to be solved analytically and hence an adaptive ODE solver COLSYS is employed. Details on COLSYS can be referred to in [37, 38]. After obtaining element dynamic stiffness matrix [ki]e, global dynamic stiffness matrixK is assembled in a direct way along the meridian with regular element location vectors.


## 4. Wittrick-Williams Algorithm

The Wittrick-Williams algorithm is a counting algorithm which gives the number of frequencies below a trial valueω∗ as(17)J=J0+JKwhereJ is the total number of frequencies below ω∗; JK=sK is the sign count of the global dynamic stiffness matrixK and equals to the number of negative elements on the diagonal of an upper triangular matrix obtained fromK by applying standard Gaussian elimination without row exchange;J0 is the number of clamped-end frequencies exceeded by ω∗ and can be obtained by accumulating J0e of each element over the whole meridian as(18)J0=∑J0eTaking advantage of the self-adaptability of COLSYS, J0e can be calculated using a substructuring method. When computing [ki]e of element (e), submesh (termed as (e^)) is defined by COLSYS as is schematically shown in Figure 4.Figure 4
Submesh by COLSYS on element (e).Applying the Wittrick-Williams algorithm again on submesh(e^), (19)J0e=∑J0e^+sKe^since COLSYS is capable of controlling error tolerance due to its self-adaptability, the number of clamped-end frequencies below ω∗ on submesh (e^) has to be zero, indicating J0(e^)=0. Otherwise, the clamped-end modes can be enlarged infinitely, resulting in unsatisfaction of any given error tolerance.For a submesh element(φ^i,φ^i+1), its dynamic stiffness matrix [k^i](e^) can be obtained by linearly combining ten solutions of (12) and (15). The element-end displacements and force vectors for submesh are denoted as d(e^) and F(e^), respectively, the submesh element stiffness matrix [k^i](e^) satisfies(20)k^ie^de^=Fe^.Applying the boundary conditions in (15) in turn again, (20) is expanded as(21)k^ie^Be^=Ce^where(22)Be^=dj=1e^dj=2e^dj=3e^dj=4e^dj=5e^dj=6e^dj=7e^dj=8e^dj=9e^dj=10e^and(23)Ce^=Fj=1e^Fj=2e^Fj=3e^Fj=4e^Fj=5e^Fj=6e^Fj=7e^Fj=8e^Fj=9e^Fj=10e^Since J0(e^)=0, [B](e^) is nonsingular (otherwise ω∗ is one of a clamped-end frequency, which is conflicted with J0(e^)=0), the submesh dynamic stiffness matrix [k^i](e^) is formulated as(24)k^ie^=Ce^Be^-1.Similarly, K(e^) can be assembled regularly from [k^i](e^), and J0e=s{K(e^)} according to (19).By employing the Wittrick-Williams algorithm repeatedly, an approximate upper and lower bound of each frequency can be predicted. Simple bisection method can then be employed to obtain exact natural frequencies. Based on the proposed theory, a FORTRAN code was written by the authors to determine natural frequencies of moderately thick shells of revolution with various meridian shapes and different boundary conditions.


## 5. Numerical Examples

Numerical examples on natural frequencies of moderately thick circular cylindrical shells, elliptical shells, and parabolic and hyperbolic shells are presented in this section using the proposed approach. Through comparison with results from existing literature, the exact dynamic stiffness formulation is proven to be reliable and accurate.


### 5.1. Circular Cylindrical Shells

For circular cylindrical shells,r is constant, φ=90∘ and R=∞, the meridian coordinate is changed to the distance coordinates along the meridian. The equation of motion in (12) is still applicable with the form of(25)Jdzds=Sz.


#### 5.1.1. Convergence ofJ0

A moderately thick circular cylindrical shell with lengthL= 20m, radiusr= 2m, thicknessh = 0.2m, Young’s modulusE=200GPa, Poisson’s ratioν=0.3, and densityρ=7800kg/m3 is examined. The shell is divided equally along the meridian.J0,JK, andJ in the W-W algorithm with different number of shell elementsne atn=1, ω∗=3000 rad/s are tabulated in Table 1. It can be observed that the solution toJ0 proposed in this paper converges to zero with the increase ofne, and the proposed submesh-based solution toJ0 presented in Section 4 is reliable.Table 1
Frequency counting of a moderately thick circular cylindrical shell with differentneunder n=1, ω∗=3000 rad/s.
ne J 0 J K J = J 0 +J K 2 16 4 20 3 15 5 20 4 12 8 20 5 10 10 20 8 8 12 20 10 10 10 20 11 0 20 20 12 0 20 20 15 0 20 20


#### 5.1.2. A Simply Supported Circular Cylindrical Shell

Khalili et al. [33] calculated natural frequencies of a thick simply supported circular cylindrical shell using a three-dimensional refined higher-order theory. The geometry of the shell isL/r=1.0 and Poisson’s ratioν=0.3. Table 2 lists the lowest three natural frequencies under differenth/r andn from both the present dynamic stiffness formulation and [33], where ω0 is a frequency parameter as ω0=(1/L)E/ρ(1+ν).m denotes the order under the correspondingn. Define the relative errorη as(26)η=Ωcurrent-ΩliteratureΩliterature×100%where Ωcurrent and Ωliterature are the frequencies obtained from the current formulation and existing literature, respectively. The relative errorη between the present formulation and [33] are also tabulated in Table 2. Sign ‘-’ denotes the obtained value from the dynamic stiffness formulation is smaller than that from literature.Table 2
Natural frequencies of a circular cylindrical shell (×ω0).
h/r m  n 1 2 3 4 5 0.1 1 Present 1.06348 0.88147 0.80659 0.89439 1.11575 Ref. [33] 1.06142 0.88158 0.80860 0.89819 1.12156 η 0.19% -0.01% -0.25% -0.42% -0.52% 2 Present 2.37850 2.72107 3.15851 3.66689 4.22897 Ref. [33] 2.37417 2.71556 3.15287 3.66166 4.22432 η 0.18% 0.20% 0.18% 0.14% 0.11% 3 Present 3.96689 4.49462 5.24988 6.14424 7.12431 Ref. [33] 3.96347 4.48736 5.23619 6.12162 7.08999 η 0.09% 0.16% 0.26% 0.37% 0.48%  0.2 1 Present 1.19118 1.09378 1.18201 1.46592 1.88124 Ref. [33] 1.18622 1.09885 1.19644 1.48941 1.91475 η 0.42% -0.46% -1.21% -1.58% -1.75% 2 Present 2.39136 2.73821 3.17697 3.68507 4.24608 Ref. [33] 2.37498 2.71703 3.15526 3.66486 4.22776 η 0.69% 0.78% 0.69% 0.55% 0.43% 3 Present 3.96726 4.49576 5.25153 6.14593 6.93506 Ref. [33] 3.95328 4.46562 5.19327 6.04694 6.97100 η 0.35% 0.67% 1.12% 1.64% -0.52%Referring to Table2, frequencies from the present dynamic stiffness formulation agree excellently well with those in [33] whenh/r=0.1, where relative errors are well less than 1.0%. Forh/r=0.2, good agreement is still achieved except some individual data points whose errors are between 1.0% and 2.0%, which is still acceptable. In general, very good agreement is achieved between the present method and [33], demonstrating the proposed dynamic stiffness formulation is exact.


### 5.2. Elliptical Shells

The meridian equation of an elliptical shell inOXZ plane is given by(27)x2a2+z2b2=1wherea andb are constants. Radius of meridian curvatureR and radius of parallel circlesr can be written with respect toφ as(28)R=a2b21+cot2⁡φa2+b2cot2⁡φ3/2(29)r=a1+b2/a2cot2⁡φ.In (28) and (29),φ is the only independent variable,R andr can be defined continuously and smoothly within each shell element along the meridian, guaranteeing an accurate solution to free vibration of elliptical shells should the proposed dynamic stiffness formulation is used.


#### 5.2.1. Clamped-Free Spherical Shell Caps

Shoulda=b in (27), an elliptical shell becomes a standard spherical shell whereR is constant andr=Rsinφ. Gautham and Ganesan [34] analysed the free vibration of moderately thick hemispherical shell caps with 30∘ and 60∘ cutout shown in Figure 5. The boundary condition is clamped-free.Figure 5
Diagram of hemispherical caps with cutout.
(a) (b)The present authors calculated the lowest three nondimensional frequency parametersΩ (Ω=ωaρ/E) under different circumferential wave numbern whenh/a=0.2 using the proposed formulation. In this example, Poisson’s ratio is 0.3. Results from both the dynamic stiffness approach and [34] are given opposite in Table 3 for comparison. It can be observed that both data agree very well and most relative errorsη are less than 1.0%. In general, most frequencies are close to results in existing literature at a 10-2 magnitude.Table 3
Ω for hemispherical caps with cutout.
n m 3 0 ∘ 6 0 ∘ Present Ref. [34] η Present Ref. [34] η 0 1 1.138 1.133 0.44% 1.276 1.289 -1.01% 2 1.309 1.297 0.93% 2.033 2.013 0.99% 3 1.420 1.427 -0.49% 2.729 2.747 -0.66%  1 1 0.751 0.753 -0.27% 1.143 1.152 -0.78% 2 1.291 1.295 -0.31% 2.348 2.337 0.47% 3 1.980 1.982 -0.10% 2.796 2.810 -0.50%  2 1 0.648 0.653 -0.77% 1.043 1.053 -0.95% 2 1.668 1.688 -1.18% 2.794 2.822 -0.99% 3 2.020 2.022 -0.10% 3.167 3.150 0.54%  3 1 1.137 1.157 -1.73% 1.141 1.154 -1.13% 2 2.062 2.087 -1.20% 3.055 3.099 -1.42% 3 2.437 2.435 0.08% 3.635 3.634 0.03%  4 1 1.682 1.719 -2.15% 1.425 1.441 -1.11% 2 2.525 2.557 -1.25% 3.356 3.407 -1.50% 3 3.115 3.107 0.26% 3.784 3.782 0.05%  5 1 2.187 2.235 -2.15% 1.832 1.853 -1.13% 2 3.103 3.151 -1.52% 3.721 3.772 -1.35% 3 3.857 3.860 -0.08% 4.020 4.025 -0.12%


#### 5.2.2. Free-Free Elliptical Shells

Shim and Kang [35] studied free-free elliptical shells (see Figure 6) with three-dimensional elastic theory and Ritz method. In this case,d=0.2a and Poisson’s ratioν=0.3. Table 4 gives the lowest three nonzero frequency parameters Ω (Ω=ωaρ/G andG is the shear modulus) of both prolate and oblate elliptical shells withn from 0 to 6 whenh/a=0.1.Table 4
Frequencies of free-free ellipsoidal shells inΩ with different b/a.
n m b/a=1/2 b/a=2 Present Ref. [35] Present Ref. [35] 0 1 1.069 1.060 1.557 1.559 2 1.658 1.620 1.652 1.652 3 2.921 2.876 1.891 1.890  1 1 1.240 1.230 1.356 1.358 2 1.977 1.933 1.620 1.616 3 2.555 2.539 1.942 1.931  2 1 0.2321 0.2223 0.1521 0.1505 2 1.445 1.419 0.9069 0.8924 3 1.767 1.771 1.443 1.423  3 1 0.5702 0.5456 0.4085 0.4066 2 2.049 2.020 1.070 1.061 3 3.121 3.108 1.654 1.639  4 1 0.9882 0.9476 0.7527 0.7509 2 2.651 2.611 1.343 1.328 3 4.080 4.049 1.961 1.936  5 1 1.473 1.416 1.175 1.175 2 3.306 3.254 1.764 1.746 3 4.948 4.912 2.417 2.385  6 1 2.017 1.944 1.670 1.672 2 4.017 3.953 2.290 2.270 3 5.819 5.787 2.980 2.946Figure 6
The diagram of a free-free elliptical shell.Referring to Table4, results from the dynamic stiffness formulation agree very well with those in [35]. This verifies that the dynamic stiffness formulation proposed in this paper is effective and accurate.Further error analysis on frequency spectrum with different value ofb/a is presented in Figure 7. It is understood that the differences between frequencies from the present method and [35] are extremely minor regardless of circumferential wave numbern. Frequency curves from the dynamic stiffness formulation and [35] are almost identical, no matter the shape of the elliptical shell is prolate or oblate.Figure 7
Comparison on frequency spectrum of elliptical shells: (a)b/a=1/2; (b)b/a=2.
(a) (b)


### 5.3. Parabolic Shells

For parabolic shells, the meridian equation inOXZ plane is(30)z=c1-xa2wherea andc are nonzero constants. Similar to the elliptical shells, radius of meridian curvatureR and the parallel circler corresponding to coordinateφ can also be defined and can be referred to mathematical handbooks.


#### 5.3.1. A Clamped-Free Parabolic Dome

Tornabene and Viola [18] calculated natural frequencies of a clamped-free parabolic dome using the Generalized Differential Quadrature (GDQ) method. The geometry and material properties of the examined parabolic dome area=4m,c=2.133m, heightH=2m, thicknessh=0.1m, densityρ=7800kg/m3, Young’s modulusE=2.1×1011N/m2, and Poisson’s ratioν=0.3. The lowest five natural frequencies calculated from the present dynamic stiffness formulation and [18] are tabulated in Table 5. It is observed that frequencies obtained from the present dynamic stiffness formulation agree very well with data from [18], showing good reliability and accuracy of the proposed formulation.Table 5
Comparison on frequencies of a clamped-free parabolic dome (Hz).
Total order n Present Ref. [18] η 1 2 91.68 91.24 0.48% 2 3 132.33 131.29 0.79% 3 4 154.79 154.91 -0.08% 4 1 160.63 157.99 1.67% 5 5 168.19 169.46 -0.75%


#### 5.3.2. Apex-Closed Parabolic Caps

Al-Khatib and Buchanan [20] examined the free vibration of clamped parabolic shell caps withν=0.3,h=1/20c. Apexes of parabolic caps are closed. A tiny circular hole with the size of φ0=0.005∘ is employed at the apexes so that singularity of differential matrix [L] in (1) can be avoided whilst not reducing accuracy. Table 6 lists the first natural frequencies of eachn whena=c. Table 7 shows the first three frequencies ofn=0 and 1 whena=2c. Frequencies are presented in the form of nondimensional parameter Ω, where Ω=ωcρ/G.Table 6
Frequencies inΩ for an apex-closed parabolic cap with a=c.
n Present Ref [20] η 0 1.798 1.774 1.35% 1 1.162 1.159 0.26% 2 1.180 1.166 1.20% 3 1.157 1.133 2.12% 4 1.283 1.240 3.47%Table 7
Frequencies inΩ for an apex-closed parabolic cap with a=2c.
n m Present Ref [20] η 0 1 0.6998 0.6955 0.62% 2 0.8004 0.7896 1.37% 3 0.9263 0.9177 0.94%  1 1 0.6058 0.6043 0.25% 2 0.7325 0.7252 1.01% 3 0.8793 0.8612 2.10%According to Tables6 and 7, good agreement is achieved between the present method and [20], demonstrating that the proposed dynamic stiffness formulation is capable of analysing the free vibration of apex-closed parabolic shells of revolution.


### 5.4. A Hyperbolic Cooling Tower Shell

Cooling towers are commonly found in civil engineering infrastructures, e.g., power plants. Their meridian can be approximated as a hyperbolic curve. Based on the illustration of Figure8, the meridian equation of a cooling tower shell inOXZ plane is given by(31)x2a2-z2b2=1.Accordingly, the radius of meridian curvatureR and parallel circler corresponding to coordinateφ are defined as(32)R=a2b21+cot2⁡φa2-b2cot2⁡φ3/2(33)r=a1-b2/a2cot2⁡φ.In this section, natural frequencies of a cooling tower at Didcot Power Station [36] are calculated. Actually the cooling tower is assumed to be a general shell of revolution; however, its basic shape can be approximated as a hyperboloid. Parameters of the hyperboloid area=25.1205m,b=60.2019m,tz=20.7264m, andbz=93.5736m. Thickness of the hyperbolic shell ish=0.1777m. The properties of material areE=2.7579×1010N/m2,ρ=2323kg/m3,ν=0.2. The boundary condition is clamped-free.Figure 8
A hyperbolic shell and its coordinate system.Strictly speaking, the cooling tower should be considered as a thin shell sinceh/a orh/b is well smaller than 1/20, which is the boundary between thin and thick shells. However, since the proposed dynamic stiffness formulation is exact, it is capable of producing accurate natural frequencies even if the shell is very thin. For comparison, results obtained from the dynamic stiffness formulation using current moderately thick shell theory, thin shell theory (previously proposed by present authors in [29]), as well as those in [36] by Deb Nath using ring shell elements are tabulated in Table 8.Table 8
Frequencies of Didcot cooling tower (Hz).
n m Thick Thin Deb Nath [36] η thick- [ 38 ] η thin- [ 38 ] 0 1 6.7463 6.7463 6.7493 -0.04% -0.04% 2 8.1912 8.1913 8.1956 -0.05% -0.05% 3 12.110 12.111 12.129 -0.16% -0.15%  1 1 3.4558 3.4560 3.4575 -0.05% -0.04% 2 7.0329 7.0331 7.0364 -0.05% -0.05% 3 11.004 11.005 11.011 -0.06% -0.05%  2 1 1.9364 1.9394 1.9373 -0.05% 0.11% 2 3.6550 3.6567 3.6571 -0.06% -0.01% 3 6.9282 6.9292 6.9320 -0.05% -0.04%  3 1 1.5935 1.6027 1.5942 -0.04% 0.53% 2 2.0114 2.0192 2.0125 -0.05% 0.33% 3 4.1627 4.1661 4.1653 -0.06% 0.02%  4 1 1.2721 1.2909 1.2730 -0.07% 1.41% 2 1.8300 1.8488 1.8308 -0.04% 0.98% 3 2.6961 2.7042 2.6979 -0.07% 0.23%  5 1 1.3206 1.3468 1.3218 -0.09% 1.89% 2 1.6820 1.7104 1.6837 -0.10% 1.59% 3 2.4524 2.4710 2.4536 -0.05% 0.71%  6 1 1.5014 1.5253 1.5030 -0.11% 1.48% 2 1.8584 1.8942 1.8604 -0.11% 1.82% 3 2.3832 2.4112 2.3861 -0.12% 1.05%  7 1 1.6851 1.7048 1.6868 -0.10% 1.07% 2 2.2235 2.2546 2.2263 -0.13% 1.27% 3 2.5470 2.5861 2.5501 -0.12% 1.41%  8 1 1.8771 1.8953 1.8792 -0.11% 0.86% 2 2.5006 2.5263 2.5038 -0.13% 0.90% 3 3.1195 3.1548 3.1236 -0.13% 1.00%Referring to Table8, conclusion can be reached that frequencies predicted both from the moderately thick and from thin dynamic stiffness formulae agree very well with data in [36]. A further comparison on the error of ηthick-[38] and ηthin-[38] shows that results obtained with moderately thick dynamic stiffness formulation enjoys a much better precision over their thin counterpart, especially when circumferential wave numbern increases. This observation suggests that although the proposed dynamic stiffness formulation is based on moderately thick shell theory, the present method can be applied to predict natural frequencies of thin shells, without observation of any shear locking. As a result, the exactness of the proposed method is validated.The actual Didcot tower is slightly different from the approximation of a hyperboloid. Figure9 plots the frequency curves calculated both as a moderately thick hyperbolic shell with present formulation and a general shell of revolution (with slightly different meridian shape based on the actual construction drawings as is explained in [36]) using FEM in [36]. It can be observed that results agree very well for differentn andm combinations, and plots for strict hyperboloid and general shell are almost identical. This further demonstrates that the assumption of a cooling tower as a hyperboloid is acceptable, and results from the proposed dynamic stiffness formulation are highly reliable.Figure 9
Comparison on frequencies obtained as a hyperboloid and a general shell.


## 5.1. Circular Cylindrical Shells

For circular cylindrical shells,r is constant, φ=90∘ and R=∞, the meridian coordinate is changed to the distance coordinates along the meridian. The equation of motion in (12) is still applicable with the form of(25)Jdzds=Sz.


### 5.1.1. Convergence ofJ0

A moderately thick circular cylindrical shell with lengthL= 20m, radiusr= 2m, thicknessh = 0.2m, Young’s modulusE=200GPa, Poisson’s ratioν=0.3, and densityρ=7800kg/m3 is examined. The shell is divided equally along the meridian.J0,JK, andJ in the W-W algorithm with different number of shell elementsne atn=1, ω∗=3000 rad/s are tabulated in Table 1. It can be observed that the solution toJ0 proposed in this paper converges to zero with the increase ofne, and the proposed submesh-based solution toJ0 presented in Section 4 is reliable.Table 1
Frequency counting of a moderately thick circular cylindrical shell with differentneunder n=1, ω∗=3000 rad/s.
ne J 0 J K J = J 0 +J K 2 16 4 20 3 15 5 20 4 12 8 20 5 10 10 20 8 8 12 20 10 10 10 20 11 0 20 20 12 0 20 20 15 0 20 20


### 5.1.2. A Simply Supported Circular Cylindrical Shell

Khalili et al. [33] calculated natural frequencies of a thick simply supported circular cylindrical shell using a three-dimensional refined higher-order theory. The geometry of the shell isL/r=1.0 and Poisson’s ratioν=0.3. Table 2 lists the lowest three natural frequencies under differenth/r andn from both the present dynamic stiffness formulation and [33], where ω0 is a frequency parameter as ω0=(1/L)E/ρ(1+ν).m denotes the order under the correspondingn. Define the relative errorη as(26)η=Ωcurrent-ΩliteratureΩliterature×100%where Ωcurrent and Ωliterature are the frequencies obtained from the current formulation and existing literature, respectively. The relative errorη between the present formulation and [33] are also tabulated in Table 2. Sign ‘-’ denotes the obtained value from the dynamic stiffness formulation is smaller than that from literature.Table 2
Natural frequencies of a circular cylindrical shell (×ω0).
h/r m  n 1 2 3 4 5 0.1 1 Present 1.06348 0.88147 0.80659 0.89439 1.11575 Ref. [33] 1.06142 0.88158 0.80860 0.89819 1.12156 η 0.19% -0.01% -0.25% -0.42% -0.52% 2 Present 2.37850 2.72107 3.15851 3.66689 4.22897 Ref. [33] 2.37417 2.71556 3.15287 3.66166 4.22432 η 0.18% 0.20% 0.18% 0.14% 0.11% 3 Present 3.96689 4.49462 5.24988 6.14424 7.12431 Ref. [33] 3.96347 4.48736 5.23619 6.12162 7.08999 η 0.09% 0.16% 0.26% 0.37% 0.48%  0.2 1 Present 1.19118 1.09378 1.18201 1.46592 1.88124 Ref. [33] 1.18622 1.09885 1.19644 1.48941 1.91475 η 0.42% -0.46% -1.21% -1.58% -1.75% 2 Present 2.39136 2.73821 3.17697 3.68507 4.24608 Ref. [33] 2.37498 2.71703 3.15526 3.66486 4.22776 η 0.69% 0.78% 0.69% 0.55% 0.43% 3 Present 3.96726 4.49576 5.25153 6.14593 6.93506 Ref. [33] 3.95328 4.46562 5.19327 6.04694 6.97100 η 0.35% 0.67% 1.12% 1.64% -0.52%Referring to Table2, frequencies from the present dynamic stiffness formulation agree excellently well with those in [33] whenh/r=0.1, where relative errors are well less than 1.0%. Forh/r=0.2, good agreement is still achieved except some individual data points whose errors are between 1.0% and 2.0%, which is still acceptable. In general, very good agreement is achieved between the present method and [33], demonstrating the proposed dynamic stiffness formulation is exact.


## 5.1.1. Convergence ofJ0

A moderately thick circular cylindrical shell with lengthL= 20m, radiusr= 2m, thicknessh = 0.2m, Young’s modulusE=200GPa, Poisson’s ratioν=0.3, and densityρ=7800kg/m3 is examined. The shell is divided equally along the meridian.J0,JK, andJ in the W-W algorithm with different number of shell elementsne atn=1, ω∗=3000 rad/s are tabulated in Table 1. It can be observed that the solution toJ0 proposed in this paper converges to zero with the increase ofne, and the proposed submesh-based solution toJ0 presented in Section 4 is reliable.Table 1
Frequency counting of a moderately thick circular cylindrical shell with differentneunder n=1, ω∗=3000 rad/s.
ne J 0 J K J = J 0 +J K 2 16 4 20 3 15 5 20 4 12 8 20 5 10 10 20 8 8 12 20 10 10 10 20 11 0 20 20 12 0 20 20 15 0 20 20


## 5.1.2. A Simply Supported Circular Cylindrical Shell

Khalili et al. [33] calculated natural frequencies of a thick simply supported circular cylindrical shell using a three-dimensional refined higher-order theory. The geometry of the shell isL/r=1.0 and Poisson’s ratioν=0.3. Table 2 lists the lowest three natural frequencies under differenth/r andn from both the present dynamic stiffness formulation and [33], where ω0 is a frequency parameter as ω0=(1/L)E/ρ(1+ν).m denotes the order under the correspondingn. Define the relative errorη as(26)η=Ωcurrent-ΩliteratureΩliterature×100%where Ωcurrent and Ωliterature are the frequencies obtained from the current formulation and existing literature, respectively. The relative errorη between the present formulation and [33] are also tabulated in Table 2. Sign ‘-’ denotes the obtained value from the dynamic stiffness formulation is smaller than that from literature.Table 2
Natural frequencies of a circular cylindrical shell (×ω0).
h/r m  n 1 2 3 4 5 0.1 1 Present 1.06348 0.88147 0.80659 0.89439 1.11575 Ref. [33] 1.06142 0.88158 0.80860 0.89819 1.12156 η 0.19% -0.01% -0.25% -0.42% -0.52% 2 Present 2.37850 2.72107 3.15851 3.66689 4.22897 Ref. [33] 2.37417 2.71556 3.15287 3.66166 4.22432 η 0.18% 0.20% 0.18% 0.14% 0.11% 3 Present 3.96689 4.49462 5.24988 6.14424 7.12431 Ref. [33] 3.96347 4.48736 5.23619 6.12162 7.08999 η 0.09% 0.16% 0.26% 0.37% 0.48%  0.2 1 Present 1.19118 1.09378 1.18201 1.46592 1.88124 Ref. [33] 1.18622 1.09885 1.19644 1.48941 1.91475 η 0.42% -0.46% -1.21% -1.58% -1.75% 2 Present 2.39136 2.73821 3.17697 3.68507 4.24608 Ref. [33] 2.37498 2.71703 3.15526 3.66486 4.22776 η 0.69% 0.78% 0.69% 0.55% 0.43% 3 Present 3.96726 4.49576 5.25153 6.14593 6.93506 Ref. [33] 3.95328 4.46562 5.19327 6.04694 6.97100 η 0.35% 0.67% 1.12% 1.64% -0.52%Referring to Table2, frequencies from the present dynamic stiffness formulation agree excellently well with those in [33] whenh/r=0.1, where relative errors are well less than 1.0%. Forh/r=0.2, good agreement is still achieved except some individual data points whose errors are between 1.0% and 2.0%, which is still acceptable. In general, very good agreement is achieved between the present method and [33], demonstrating the proposed dynamic stiffness formulation is exact.


## 5.2. Elliptical Shells

The meridian equation of an elliptical shell inOXZ plane is given by(27)x2a2+z2b2=1wherea andb are constants. Radius of meridian curvatureR and radius of parallel circlesr can be written with respect toφ as(28)R=a2b21+cot2⁡φa2+b2cot2⁡φ3/2(29)r=a1+b2/a2cot2⁡φ.In (28) and (29),φ is the only independent variable,R andr can be defined continuously and smoothly within each shell element along the meridian, guaranteeing an accurate solution to free vibration of elliptical shells should the proposed dynamic stiffness formulation is used.


### 5.2.1. Clamped-Free Spherical Shell Caps

Shoulda=b in (27), an elliptical shell becomes a standard spherical shell whereR is constant andr=Rsinφ. Gautham and Ganesan [34] analysed the free vibration of moderately thick hemispherical shell caps with 30∘ and 60∘ cutout shown in Figure 5. The boundary condition is clamped-free.Figure 5
Diagram of hemispherical caps with cutout.
(a) (b)The present authors calculated the lowest three nondimensional frequency parametersΩ (Ω=ωaρ/E) under different circumferential wave numbern whenh/a=0.2 using the proposed formulation. In this example, Poisson’s ratio is 0.3. Results from both the dynamic stiffness approach and [34] are given opposite in Table 3 for comparison. It can be observed that both data agree very well and most relative errorsη are less than 1.0%. In general, most frequencies are close to results in existing literature at a 10-2 magnitude.Table 3
Ω for hemispherical caps with cutout.
n m 3 0 ∘ 6 0 ∘ Present Ref. [34] η Present Ref. [34] η 0 1 1.138 1.133 0.44% 1.276 1.289 -1.01% 2 1.309 1.297 0.93% 2.033 2.013 0.99% 3 1.420 1.427 -0.49% 2.729 2.747 -0.66%  1 1 0.751 0.753 -0.27% 1.143 1.152 -0.78% 2 1.291 1.295 -0.31% 2.348 2.337 0.47% 3 1.980 1.982 -0.10% 2.796 2.810 -0.50%  2 1 0.648 0.653 -0.77% 1.043 1.053 -0.95% 2 1.668 1.688 -1.18% 2.794 2.822 -0.99% 3 2.020 2.022 -0.10% 3.167 3.150 0.54%  3 1 1.137 1.157 -1.73% 1.141 1.154 -1.13% 2 2.062 2.087 -1.20% 3.055 3.099 -1.42% 3 2.437 2.435 0.08% 3.635 3.634 0.03%  4 1 1.682 1.719 -2.15% 1.425 1.441 -1.11% 2 2.525 2.557 -1.25% 3.356 3.407 -1.50% 3 3.115 3.107 0.26% 3.784 3.782 0.05%  5 1 2.187 2.235 -2.15% 1.832 1.853 -1.13% 2 3.103 3.151 -1.52% 3.721 3.772 -1.35% 3 3.857 3.860 -0.08% 4.020 4.025 -0.12%


### 5.2.2. Free-Free Elliptical Shells

Shim and Kang [35] studied free-free elliptical shells (see Figure 6) with three-dimensional elastic theory and Ritz method. In this case,d=0.2a and Poisson’s ratioν=0.3. Table 4 gives the lowest three nonzero frequency parameters Ω (Ω=ωaρ/G andG is the shear modulus) of both prolate and oblate elliptical shells withn from 0 to 6 whenh/a=0.1.Table 4
Frequencies of free-free ellipsoidal shells inΩ with different b/a.
n m b/a=1/2 b/a=2 Present Ref. [35] Present Ref. [35] 0 1 1.069 1.060 1.557 1.559 2 1.658 1.620 1.652 1.652 3 2.921 2.876 1.891 1.890  1 1 1.240 1.230 1.356 1.358 2 1.977 1.933 1.620 1.616 3 2.555 2.539 1.942 1.931  2 1 0.2321 0.2223 0.1521 0.1505 2 1.445 1.419 0.9069 0.8924 3 1.767 1.771 1.443 1.423  3 1 0.5702 0.5456 0.4085 0.4066 2 2.049 2.020 1.070 1.061 3 3.121 3.108 1.654 1.639  4 1 0.9882 0.9476 0.7527 0.7509 2 2.651 2.611 1.343 1.328 3 4.080 4.049 1.961 1.936  5 1 1.473 1.416 1.175 1.175 2 3.306 3.254 1.764 1.746 3 4.948 4.912 2.417 2.385  6 1 2.017 1.944 1.670 1.672 2 4.017 3.953 2.290 2.270 3 5.819 5.787 2.980 2.946Figure 6
The diagram of a free-free elliptical shell.Referring to Table4, results from the dynamic stiffness formulation agree very well with those in [35]. This verifies that the dynamic stiffness formulation proposed in this paper is effective and accurate.Further error analysis on frequency spectrum with different value ofb/a is presented in Figure 7. It is understood that the differences between frequencies from the present method and [35] are extremely minor regardless of circumferential wave numbern. Frequency curves from the dynamic stiffness formulation and [35] are almost identical, no matter the shape of the elliptical shell is prolate or oblate.Figure 7
Comparison on frequency spectrum of elliptical shells: (a)b/a=1/2; (b)b/a=2.
(a) (b)


## 5.2.1. Clamped-Free Spherical Shell Caps

Shoulda=b in (27), an elliptical shell becomes a standard spherical shell whereR is constant andr=Rsinφ. Gautham and Ganesan [34] analysed the free vibration of moderately thick hemispherical shell caps with 30∘ and 60∘ cutout shown in Figure 5. The boundary condition is clamped-free.Figure 5
Diagram of hemispherical caps with cutout.
(a) (b)The present authors calculated the lowest three nondimensional frequency parametersΩ (Ω=ωaρ/E) under different circumferential wave numbern whenh/a=0.2 using the proposed formulation. In this example, Poisson’s ratio is 0.3. Results from both the dynamic stiffness approach and [34] are given opposite in Table 3 for comparison. It can be observed that both data agree very well and most relative errorsη are less than 1.0%. In general, most frequencies are close to results in existing literature at a 10-2 magnitude.Table 3
Ω for hemispherical caps with cutout.
n m 3 0 ∘ 6 0 ∘ Present Ref. [34] η Present Ref. [34] η 0 1 1.138 1.133 0.44% 1.276 1.289 -1.01% 2 1.309 1.297 0.93% 2.033 2.013 0.99% 3 1.420 1.427 -0.49% 2.729 2.747 -0.66%  1 1 0.751 0.753 -0.27% 1.143 1.152 -0.78% 2 1.291 1.295 -0.31% 2.348 2.337 0.47% 3 1.980 1.982 -0.10% 2.796 2.810 -0.50%  2 1 0.648 0.653 -0.77% 1.043 1.053 -0.95% 2 1.668 1.688 -1.18% 2.794 2.822 -0.99% 3 2.020 2.022 -0.10% 3.167 3.150 0.54%  3 1 1.137 1.157 -1.73% 1.141 1.154 -1.13% 2 2.062 2.087 -1.20% 3.055 3.099 -1.42% 3 2.437 2.435 0.08% 3.635 3.634 0.03%  4 1 1.682 1.719 -2.15% 1.425 1.441 -1.11% 2 2.525 2.557 -1.25% 3.356 3.407 -1.50% 3 3.115 3.107 0.26% 3.784 3.782 0.05%  5 1 2.187 2.235 -2.15% 1.832 1.853 -1.13% 2 3.103 3.151 -1.52% 3.721 3.772 -1.35% 3 3.857 3.860 -0.08% 4.020 4.025 -0.12%


## 5.2.2. Free-Free Elliptical Shells

Shim and Kang [35] studied free-free elliptical shells (see Figure 6) with three-dimensional elastic theory and Ritz method. In this case,d=0.2a and Poisson’s ratioν=0.3. Table 4 gives the lowest three nonzero frequency parameters Ω (Ω=ωaρ/G andG is the shear modulus) of both prolate and oblate elliptical shells withn from 0 to 6 whenh/a=0.1.Table 4
Frequencies of free-free ellipsoidal shells inΩ with different b/a.
n m b/a=1/2 b/a=2 Present Ref. [35] Present Ref. [35] 0 1 1.069 1.060 1.557 1.559 2 1.658 1.620 1.652 1.652 3 2.921 2.876 1.891 1.890  1 1 1.240 1.230 1.356 1.358 2 1.977 1.933 1.620 1.616 3 2.555 2.539 1.942 1.931  2 1 0.2321 0.2223 0.1521 0.1505 2 1.445 1.419 0.9069 0.8924 3 1.767 1.771 1.443 1.423  3 1 0.5702 0.5456 0.4085 0.4066 2 2.049 2.020 1.070 1.061 3 3.121 3.108 1.654 1.639  4 1 0.9882 0.9476 0.7527 0.7509 2 2.651 2.611 1.343 1.328 3 4.080 4.049 1.961 1.936  5 1 1.473 1.416 1.175 1.175 2 3.306 3.254 1.764 1.746 3 4.948 4.912 2.417 2.385  6 1 2.017 1.944 1.670 1.672 2 4.017 3.953 2.290 2.270 3 5.819 5.787 2.980 2.946Figure 6
The diagram of a free-free elliptical shell.Referring to Table4, results from the dynamic stiffness formulation agree very well with those in [35]. This verifies that the dynamic stiffness formulation proposed in this paper is effective and accurate.Further error analysis on frequency spectrum with different value ofb/a is presented in Figure 7. It is understood that the differences between frequencies from the present method and [35] are extremely minor regardless of circumferential wave numbern. Frequency curves from the dynamic stiffness formulation and [35] are almost identical, no matter the shape of the elliptical shell is prolate or oblate.Figure 7
Comparison on frequency spectrum of elliptical shells: (a)b/a=1/2; (b)b/a=2.
(a) (b)


## 5.3. Parabolic Shells

For parabolic shells, the meridian equation inOXZ plane is(30)z=c1-xa2wherea andc are nonzero constants. Similar to the elliptical shells, radius of meridian curvatureR and the parallel circler corresponding to coordinateφ can also be defined and can be referred to mathematical handbooks.


### 5.3.1. A Clamped-Free Parabolic Dome

Tornabene and Viola [18] calculated natural frequencies of a clamped-free parabolic dome using the Generalized Differential Quadrature (GDQ) method. The geometry and material properties of the examined parabolic dome area=4m,c=2.133m, heightH=2m, thicknessh=0.1m, densityρ=7800kg/m3, Young’s modulusE=2.1×1011N/m2, and Poisson’s ratioν=0.3. The lowest five natural frequencies calculated from the present dynamic stiffness formulation and [18] are tabulated in Table 5. It is observed that frequencies obtained from the present dynamic stiffness formulation agree very well with data from [18], showing good reliability and accuracy of the proposed formulation.Table 5
Comparison on frequencies of a clamped-free parabolic dome (Hz).
Total order n Present Ref. [18] η 1 2 91.68 91.24 0.48% 2 3 132.33 131.29 0.79% 3 4 154.79 154.91 -0.08% 4 1 160.63 157.99 1.67% 5 5 168.19 169.46 -0.75%


### 5.3.2. Apex-Closed Parabolic Caps

Al-Khatib and Buchanan [20] examined the free vibration of clamped parabolic shell caps withν=0.3,h=1/20c. Apexes of parabolic caps are closed. A tiny circular hole with the size of φ0=0.005∘ is employed at the apexes so that singularity of differential matrix [L] in (1) can be avoided whilst not reducing accuracy. Table 6 lists the first natural frequencies of eachn whena=c. Table 7 shows the first three frequencies ofn=0 and 1 whena=2c. Frequencies are presented in the form of nondimensional parameter Ω, where Ω=ωcρ/G.Table 6
Frequencies inΩ for an apex-closed parabolic cap with a=c.
n Present Ref [20] η 0 1.798 1.774 1.35% 1 1.162 1.159 0.26% 2 1.180 1.166 1.20% 3 1.157 1.133 2.12% 4 1.283 1.240 3.47%Table 7
Frequencies inΩ for an apex-closed parabolic cap with a=2c.
n m Present Ref [20] η 0 1 0.6998 0.6955 0.62% 2 0.8004 0.7896 1.37% 3 0.9263 0.9177 0.94%  1 1 0.6058 0.6043 0.25% 2 0.7325 0.7252 1.01% 3 0.8793 0.8612 2.10%According to Tables6 and 7, good agreement is achieved between the present method and [20], demonstrating that the proposed dynamic stiffness formulation is capable of analysing the free vibration of apex-closed parabolic shells of revolution.


## 5.3.1. A Clamped-Free Parabolic Dome

Tornabene and Viola [18] calculated natural frequencies of a clamped-free parabolic dome using the Generalized Differential Quadrature (GDQ) method. The geometry and material properties of the examined parabolic dome area=4m,c=2.133m, heightH=2m, thicknessh=0.1m, densityρ=7800kg/m3, Young’s modulusE=2.1×1011N/m2, and Poisson’s ratioν=0.3. The lowest five natural frequencies calculated from the present dynamic stiffness formulation and [18] are tabulated in Table 5. It is observed that frequencies obtained from the present dynamic stiffness formulation agree very well with data from [18], showing good reliability and accuracy of the proposed formulation.Table 5
Comparison on frequencies of a clamped-free parabolic dome (Hz).
Total order n Present Ref. [18] η 1 2 91.68 91.24 0.48% 2 3 132.33 131.29 0.79% 3 4 154.79 154.91 -0.08% 4 1 160.63 157.99 1.67% 5 5 168.19 169.46 -0.75%


## 5.3.2. Apex-Closed Parabolic Caps

Al-Khatib and Buchanan [20] examined the free vibration of clamped parabolic shell caps withν=0.3,h=1/20c. Apexes of parabolic caps are closed. A tiny circular hole with the size of φ0=0.005∘ is employed at the apexes so that singularity of differential matrix [L] in (1) can be avoided whilst not reducing accuracy. Table 6 lists the first natural frequencies of eachn whena=c. Table 7 shows the first three frequencies ofn=0 and 1 whena=2c. Frequencies are presented in the form of nondimensional parameter Ω, where Ω=ωcρ/G.Table 6
Frequencies inΩ for an apex-closed parabolic cap with a=c.
n Present Ref [20] η 0 1.798 1.774 1.35% 1 1.162 1.159 0.26% 2 1.180 1.166 1.20% 3 1.157 1.133 2.12% 4 1.283 1.240 3.47%Table 7
Frequencies inΩ for an apex-closed parabolic cap with a=2c.
n m Present Ref [20] η 0 1 0.6998 0.6955 0.62% 2 0.8004 0.7896 1.37% 3 0.9263 0.9177 0.94%  1 1 0.6058 0.6043 0.25% 2 0.7325 0.7252 1.01% 3 0.8793 0.8612 2.10%According to Tables6 and 7, good agreement is achieved between the present method and [20], demonstrating that the proposed dynamic stiffness formulation is capable of analysing the free vibration of apex-closed parabolic shells of revolution.


## 5.4. A Hyperbolic Cooling Tower Shell

Cooling towers are commonly found in civil engineering infrastructures, e.g., power plants. Their meridian can be approximated as a hyperbolic curve. Based on the illustration of Figure8, the meridian equation of a cooling tower shell inOXZ plane is given by(31)x2a2-z2b2=1.Accordingly, the radius of meridian curvatureR and parallel circler corresponding to coordinateφ are defined as(32)R=a2b21+cot2⁡φa2-b2cot2⁡φ3/2(33)r=a1-b2/a2cot2⁡φ.In this section, natural frequencies of a cooling tower at Didcot Power Station [36] are calculated. Actually the cooling tower is assumed to be a general shell of revolution; however, its basic shape can be approximated as a hyperboloid. Parameters of the hyperboloid area=25.1205m,b=60.2019m,tz=20.7264m, andbz=93.5736m. Thickness of the hyperbolic shell ish=0.1777m. The properties of material areE=2.7579×1010N/m2,ρ=2323kg/m3,ν=0.2. The boundary condition is clamped-free.Figure 8
A hyperbolic shell and its coordinate system.Strictly speaking, the cooling tower should be considered as a thin shell sinceh/a orh/b is well smaller than 1/20, which is the boundary between thin and thick shells. However, since the proposed dynamic stiffness formulation is exact, it is capable of producing accurate natural frequencies even if the shell is very thin. For comparison, results obtained from the dynamic stiffness formulation using current moderately thick shell theory, thin shell theory (previously proposed by present authors in [29]), as well as those in [36] by Deb Nath using ring shell elements are tabulated in Table 8.Table 8
Frequencies of Didcot cooling tower (Hz).
n m Thick Thin Deb Nath [36] η thick- [ 38 ] η thin- [ 38 ] 0 1 6.7463 6.7463 6.7493 -0.04% -0.04% 2 8.1912 8.1913 8.1956 -0.05% -0.05% 3 12.110 12.111 12.129 -0.16% -0.15%  1 1 3.4558 3.4560 3.4575 -0.05% -0.04% 2 7.0329 7.0331 7.0364 -0.05% -0.05% 3 11.004 11.005 11.011 -0.06% -0.05%  2 1 1.9364 1.9394 1.9373 -0.05% 0.11% 2 3.6550 3.6567 3.6571 -0.06% -0.01% 3 6.9282 6.9292 6.9320 -0.05% -0.04%  3 1 1.5935 1.6027 1.5942 -0.04% 0.53% 2 2.0114 2.0192 2.0125 -0.05% 0.33% 3 4.1627 4.1661 4.1653 -0.06% 0.02%  4 1 1.2721 1.2909 1.2730 -0.07% 1.41% 2 1.8300 1.8488 1.8308 -0.04% 0.98% 3 2.6961 2.7042 2.6979 -0.07% 0.23%  5 1 1.3206 1.3468 1.3218 -0.09% 1.89% 2 1.6820 1.7104 1.6837 -0.10% 1.59% 3 2.4524 2.4710 2.4536 -0.05% 0.71%  6 1 1.5014 1.5253 1.5030 -0.11% 1.48% 2 1.8584 1.8942 1.8604 -0.11% 1.82% 3 2.3832 2.4112 2.3861 -0.12% 1.05%  7 1 1.6851 1.7048 1.6868 -0.10% 1.07% 2 2.2235 2.2546 2.2263 -0.13% 1.27% 3 2.5470 2.5861 2.5501 -0.12% 1.41%  8 1 1.8771 1.8953 1.8792 -0.11% 0.86% 2 2.5006 2.5263 2.5038 -0.13% 0.90% 3 3.1195 3.1548 3.1236 -0.13% 1.00%Referring to Table8, conclusion can be reached that frequencies predicted both from the moderately thick and from thin dynamic stiffness formulae agree very well with data in [36]. A further comparison on the error of ηthick-[38] and ηthin-[38] shows that results obtained with moderately thick dynamic stiffness formulation enjoys a much better precision over their thin counterpart, especially when circumferential wave numbern increases. This observation suggests that although the proposed dynamic stiffness formulation is based on moderately thick shell theory, the present method can be applied to predict natural frequencies of thin shells, without observation of any shear locking. As a result, the exactness of the proposed method is validated.The actual Didcot tower is slightly different from the approximation of a hyperboloid. Figure9 plots the frequency curves calculated both as a moderately thick hyperbolic shell with present formulation and a general shell of revolution (with slightly different meridian shape based on the actual construction drawings as is explained in [36]) using FEM in [36]. It can be observed that results agree very well for differentn andm combinations, and plots for strict hyperboloid and general shell are almost identical. This further demonstrates that the assumption of a cooling tower as a hyperboloid is acceptable, and results from the proposed dynamic stiffness formulation are highly reliable.Figure 9
Comparison on frequencies obtained as a hyperboloid and a general shell.


## 6. Further Discussions

In this section, a further study is conducted to investigate the influences of parameters (e.g., radius, thickness, material properties, etc.) on natural frequencies of cylindrical and spherical shells.


### 6.1. Cylindrical Shells

Consider a clamped-clamped (C-C) circular cylindrical shell with lengthL=20m. Baseline parameters are radiusr=2m, thicknessh=0.2m, Poisson’s ratioν=0.3, Young’s modulusE=200GPa, and densityρ=7800kg/m3. When one of these parameters changes, the rest of them are kept unaltered. To better illustrate the inherent relations between frequencies and shell parameters, a nondimensional frequency parameter Ω is defined as(34)Ω=ωrρ1-ν2E.Therefore, irrelevant parameters such as Young’s modulusE and densityρ to Ω can be ignored.Figure10 shows the trends of the lowest natural frequency parameter Ω with circumferential wave numbern from 0 to 4 against radiusr. Referring to Figure 10, frequencies with differentn exhibit differently with the increase of radiusr. Whenn is small, e.g.,n=0 and 1, the lowest frequency increases almost linearly with the increase ofr. However, this trend is not true any longer forn=2, 3, and 4, as the lowest natural frequency for a smallerr is much higher than that whenr increases.Figure 10
Trends of the lowest natural frequency parameter Ω against radiusr.Figure11 plots the trends of Ω against thicknessh. It can be observed that, for a smalln (n=0 and 1), increase of thicknessh has little effect on the value of Ω, and both curves are almost horizontal. However, a monotonic linear increase trend is obtained forn=2, 3, and 4. A possible conclusion can be reached that vibration modes withn=0 and 1 is not sensitive to thicknessh.Figure 11
Trends of the lowest natural frequency parameter Ω against thicknessh.Influence of Poisson’s ratioν on frequency parameter Ω is given in Figure 12. Referring to Figure 12, a minor effect from Poisson’s ratioν on the frequency parameter Ω is obtained when it varies.Figure 12
Trends of the lowest natural frequency parameter Ω against Poisson’s ratioν.


### 6.2. Spherical Shells

Here a clamped-free (C-F) spherical shell is considered (Figure13). Same as the cylindrical shell studied in the previous section, baseline geometry and material parameters for this spherical shell are defined as follows: radiusa=1.0m, thicknessh=0.15m (Figure 15), densityρ=2700kg/m3, Poisson’s ratioν=0.3, and Young’s modulusE=70GPa.Figure 13
Illustration of a clamped-free spherical shell.Engineering frequenciesf (in Hz) are computed in this section and thus both Young’s modulusE and densityρ play a role in the obtained natural frequencies. Parametersa,h,ρ, andE are investigated. Figure 14 shows the change of natural frequencyf against radiusa. The frequency decreases dramatically whena increases, and a converged value with differentn can be expected whena is large enough. For thef-h relation, it is quite similar to that in the previously examined circular cylindrical shell: whenn is small (n=0 and 1), the influence is almost neglectable; whenn is large (n=2, 3 and 4),f increases with the increase ofh.Figure 14
Trends of natural frequencies of a spherical shell against radiusa.Figure 15
Trends of natural frequencies of a spherical shell against thicknessh.Densityρ and Young’s modulusE are vital to engineering frequencyf. Figure 16 shows a decrease trend off with the increase of densityρ. This can be simply interpreted. With the increase of densityρ, the mass matrix increases accordingly. Since the stiffness matrix is not changed, frequency decreases naturally according to the general eigenvalue definitions. Regardingf-E relation, plotting is given in Figure 17 and similar conclusion applies to the spherical shell. Furthermore, it is worth mentioning that the smallest frequency of a shell is not necessarily obtained with a smallestn, i.e.,n=0. On the contrary, it is determined by the vibration mode that is the easiest to occur, which is usuallyn=2 in this example and other shell analysis with normal dimensions. It can also be summarised that the inherent characteristics of free vibration of moderately thick shells of revolution is revealed by the proposed dynamic stiffness formulation.Figure 16
Trends of natural frequencies of a spherical shell against densityρ.Figure 17
Trends of natural frequencies of a spherical shell against Young’s modulusE.


## 6.1. Cylindrical Shells

Consider a clamped-clamped (C-C) circular cylindrical shell with lengthL=20m. Baseline parameters are radiusr=2m, thicknessh=0.2m, Poisson’s ratioν=0.3, Young’s modulusE=200GPa, and densityρ=7800kg/m3. When one of these parameters changes, the rest of them are kept unaltered. To better illustrate the inherent relations between frequencies and shell parameters, a nondimensional frequency parameter Ω is defined as(34)Ω=ωrρ1-ν2E.Therefore, irrelevant parameters such as Young’s modulusE and densityρ to Ω can be ignored.Figure10 shows the trends of the lowest natural frequency parameter Ω with circumferential wave numbern from 0 to 4 against radiusr. Referring to Figure 10, frequencies with differentn exhibit differently with the increase of radiusr. Whenn is small, e.g.,n=0 and 1, the lowest frequency increases almost linearly with the increase ofr. However, this trend is not true any longer forn=2, 3, and 4, as the lowest natural frequency for a smallerr is much higher than that whenr increases.Figure 10
Trends of the lowest natural frequency parameter Ω against radiusr.Figure11 plots the trends of Ω against thicknessh. It can be observed that, for a smalln (n=0 and 1), increase of thicknessh has little effect on the value of Ω, and both curves are almost horizontal. However, a monotonic linear increase trend is obtained forn=2, 3, and 4. A possible conclusion can be reached that vibration modes withn=0 and 1 is not sensitive to thicknessh.Figure 11
Trends of the lowest natural frequency parameter Ω against thicknessh.Influence of Poisson’s ratioν on frequency parameter Ω is given in Figure 12. Referring to Figure 12, a minor effect from Poisson’s ratioν on the frequency parameter Ω is obtained when it varies.Figure 12
Trends of the lowest natural frequency parameter Ω against Poisson’s ratioν.


## 6.2. Spherical Shells

Here a clamped-free (C-F) spherical shell is considered (Figure13). Same as the cylindrical shell studied in the previous section, baseline geometry and material parameters for this spherical shell are defined as follows: radiusa=1.0m, thicknessh=0.15m (Figure 15), densityρ=2700kg/m3, Poisson’s ratioν=0.3, and Young’s modulusE=70GPa.Figure 13
Illustration of a clamped-free spherical shell.Engineering frequenciesf (in Hz) are computed in this section and thus both Young’s modulusE and densityρ play a role in the obtained natural frequencies. Parametersa,h,ρ, andE are investigated. Figure 14 shows the change of natural frequencyf against radiusa. The frequency decreases dramatically whena increases, and a converged value with differentn can be expected whena is large enough. For thef-h relation, it is quite similar to that in the previously examined circular cylindrical shell: whenn is small (n=0 and 1), the influence is almost neglectable; whenn is large (n=2, 3 and 4),f increases with the increase ofh.Figure 14
Trends of natural frequencies of a spherical shell against radiusa.Figure 15
Trends of natural frequencies of a spherical shell against thicknessh.Densityρ and Young’s modulusE are vital to engineering frequencyf. Figure 16 shows a decrease trend off with the increase of densityρ. This can be simply interpreted. With the increase of densityρ, the mass matrix increases accordingly. Since the stiffness matrix is not changed, frequency decreases naturally according to the general eigenvalue definitions. Regardingf-E relation, plotting is given in Figure 17 and similar conclusion applies to the spherical shell. Furthermore, it is worth mentioning that the smallest frequency of a shell is not necessarily obtained with a smallestn, i.e.,n=0. On the contrary, it is determined by the vibration mode that is the easiest to occur, which is usuallyn=2 in this example and other shell analysis with normal dimensions. It can also be summarised that the inherent characteristics of free vibration of moderately thick shells of revolution is revealed by the proposed dynamic stiffness formulation.Figure 16
Trends of natural frequencies of a spherical shell against densityρ.Figure 17
Trends of natural frequencies of a spherical shell against Young’s modulusE.


## 7. Conclusions

Free vibration behaviour and natural frequencies of shells of revolution are investigated using an exact dynamic stiffness formulation with the consideration of transverse shear deformation and rotary inertia. The W-W algorithm is employed onto the dynamic stiffnesses and natural frequencies are obtained accurately. Since the W-W algorithm is mathematically strict, no eigenvalue will be omitted. The number of clamped-end frequenciesJ0 in the W-W algorithm is given, proving the proposed formulation is complete and converged. When computing dynamic stiffnesses, governing vibration PDEs are degraded into a set of ODEs analytically, and a spline ODE solver COLSYS is used to calculate the dynamic stiffness matrices numerically. Therefore, the accuracy of obtained natural frequencies can be guaranteed and the proposed method is exact.Numerical examples on natural frequencies of a variety moderately thick shells of revolution including circular cylindrical, spherical, elliptical, and parabolic and hyperbolic shells are presented. Shells with different boundary conditions such as simply supported, clamped, and free ends are modelled. Results are compared with data in existing literature, and very good agreement is achieved. In general, the exact dynamic stiffness formulation for free vibration of moderately thick shells of revolution presented in this paper is reliable and accurate. Though equations are derived from Reissner-Mindlin assumptions, the formulation well applies to very thin shells. Furthermore, it outperforms its counterpart based on thin shell theory when comparing with data in open literature, without observation of shear locking.Further discussions on the influences of cylindrical and spherical shell parameters such as radius, thickness, Young’s modulus, Poisson’s ratio, and density are presented. It is found that the smallest natural frequency is determined by the vibration mode which is the easiest to occur (usually obtained withn=2 for shells with a normal dimension). The proposed dynamic stiffness formulation is proven to be a useful tool for analysing free vibration characteristics of moderately thick shells of revolution.



---

*Source: 4156360-2018-06-04.xml*