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

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

---

## Abstract

An exact generalised formulation for the free vibration of shells of revolution with general shaped meridians 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 frequenciesJ
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.


---

## Body

## 1. 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–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.


## 2. 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 Figure1,α is the angle between the normal directionn
→ of the meridian at an arbitrary point P and the positive direction of z-axis;β is the angle between the meridian plane at point P and a base meridian planeOXZ; O
1 and O
2 are the centres of the meridian curvature and the parallel circle;r and R are the radii of the parallel circle and meridian curvature at pointP, respectively. Young’s modulus of the material isE, Poisson’s ratio isν, density isρ, and the thickness ish.Figure 1
Geometry of a shell of revolution.According to the thin shell theory, the strain-displacement relationship of a shell of revolution is(1)
ε
=
L
Δ
,where the strain vectorε
=
ε
α
,
ε
β
,
ε
α
β
,
χ
α
,
χ
β
,
χ
α
β
T and the displacement vectorΔ
=
u
,
v
,
w
T. The expression of [L] can be found in [5]. Similarly, the force-strain relationship of a shell of revolution is(2)
N
=
D
ε
,where the internal force vector N
=
N
α
,
N
β
,
N
α
β
,
M
α
,
M
β
,
M
α
β
T and [
D
] is the stiffness matrix given by(3)
D
=
E
h
1
-
ν
2
1
ν
0
0
0
0
1
0
0
0
0
1
-
ν
2
0
0
0
h
2
12
0
0
s
y
m
.
h
2
12
0
h
2
1
-
ν
24
.Denote bydA 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 d
A
=
r
R
d
α
d
β. The total strain energy U for free vibration of the shell of revolution can be integrated along the whole middle surfaceS: (4)
U
=
∬
S
1
2
ε
T
D
ε
d
A
,and the total free vibration kinetic energy is obtained as(5)
T
=
ρ
2
∬
S
h
∂
u
∂
t
2
+
∂
v
∂
t
2
+
∂
w
∂
t
2
d
A
.Assume that the dynamic displacement functionsu,v, and w are(6)
u
α
,
β
,
t
v
α
,
β
,
t
w
α
,
β
,
t
=
u
α
cos
⁡
n
β
sin
⁡
ω
t
v
α
sin
⁡
n
β
sin
⁡
ω
t
w
α
cos
⁡
n
β
sin
⁡
ω
t
n
=
0,1
,
2
,
…
,where n is the circumferential wave number and ω is a circular frequency. According to Hamilton’s principle, (7)
δ
∫
t
1
t
2
U
-
T
d
t
=
0
.Since time interval (t
1, t
2) is arbitrary, the symplectic components (generalised displacements and forces) can be obtained from (7) after displacement substitution and integration manipulation.Figure 2
An infinitesimal portion of a shell of revolution.Referring to Figure3(a), when integrating along meridian s
→ of a shell of revolution with a positive Gaussian curvature, α
1
<
α
2, R
>
0, and r
>
0. The generalised displacements and forces obtained from (7) are given by(8)
u
1
=
u
n
,
u
2
=
v
n
,
u
3
=
w
n
,
u
4
=
w
n
′
R
,
v
1
=
r
N
-
α
,
v
2
=
r
N
-
α
β
,
v
3
=
1
R
r
M
-
α
′
-
M
-
β
cos
⁡
α
+
2
M
-
α
β
n
,
v
4
=
-
r
M
-
α
,where(
)
′ denotes differentiation with respect to α and(9)
N
-
α
=
K
ν
cos
⁡
α
r
u
1
+
u
1
′
R
+
n
ν
r
u
2
+
1
R
+
ν
sin
⁡
α
r
u
3
,
N
-
β
=
K
cos
⁡
α
r
u
1
+
ν
R
u
1
′
+
n
r
u
2
+
ν
R
+
sin
⁡
α
r
u
3
,
N
-
α
β
=
K
1
-
ν
2
-
n
r
u
1
+
r
R
u
2
r
′
,
M
-
α
=
D
ν
n
2
r
2
u
3
-
ν
cos
⁡
α
r
R
u
3
′
-
1
R
u
3
′
R
′
,
M
-
β
=
D
n
2
r
2
u
3
-
cos
⁡
α
r
R
u
3
′
-
ν
R
u
3
′
R
′
,
M
-
α
β
=
D
1
-
ν
n
r
-
cos
⁡
α
r
u
3
+
u
3
′
R
.The governing equations of motion in free vibration in the Hamilton form are given by(10)
J
z
′
=
S
z
,where z
=
u
1
u
2
u
3
u
4
v
1
v
2
v
3
v
4
T, J
=
0
-
I
I
0,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, R
<
0, and r
>
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.Figure 3
Change of the angleα according to different Gaussian curvatures.
(a)
Positive Gaussian curvature
(b)
Negative Gaussian curvature


## 3. Dynamic Stiffness Matrices

Should the circumferential wave numbern 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 (u
1
,
u
2
,
u
3
,
u
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 Figure4, 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α
aandα
b, and the segment-end displacement vector{
d
}
eand force vector{
F
}
eare defined as(11)
d
e
=
u
1
α
a
,
u
2
α
a
,
u
3
α
a
,
u
4
α
a
,
u
1
α
b
,
u
2
α
b
,
u
3
α
b
,
u
4
α
b
T
,
F
e
=
-
v
1
α
a
,
-
v
2
α
a
,
-
v
3
α
a
,
-
v
4
α
a
,
v
1
α
b
,
v
2
α
b
,
v
3
α
b
,
v
4
α
b
T
.Figure 4
Mesh divisions along the meridian.Mathematically apply the boundary conditions in (12) to (10) in turns: (12)
d
e
=
e
j
,
j
=
1
,
…
,
8
,where e
j is a unit vector with thejth element equal to one, and the corresponding element dynamic stiffness matrix [
k
i
]
e is formulated by(13)
k
i
e
=
F
j
=
1
e
F
j
=
2
e
F
j
=
3
e
F
j
=
4
e
F
j
=
5
e
F
j
=
6
e
F
j
=
7
e
F
j
=
8
e
.Due to the complexity, (10) is solved numerically using adaptive ordinary differential equations (ODEs) solver COLSYS [33, 34]. The global dynamic stiffness matrixK of a shell of revolution is assembled in a direct way with regular element location vectors.


## 4. 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(14)
J
=
J
0
+
J
K
,where J is the total number of frequencies exceeded by a specific trial frequency ω
∗, J
K
=
s
K
ω
∗ is the sign count of the global dynamic stiffness matrixK and equals the number of negative elements on the diagonal of an upper triangular matrix obtained fromK by applying standard Gaussian elimination without row exchange, and J
0 is the number of clamped-end frequencies exceeded by ω
∗ and can be accumulated from J
0
e of each shell segment element over the whole meridian as is expressed in(15)
J
0
=
∑
J
0
e
.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  (
e
^
)) are defined by COLSYS as is shown in Figure 5.Figure 5
Submesh divisions of COLSYS on element (e).Applying the W-W algorithm again on the submesh(
e
^
),(16)
J
0
e
=
∑
J
0
e
^
+
s
K
e
^
ω
∗
.Since COLSYS is capable of controlling the error tolerance due to its self-adaptability, the number of clamped-end frequencies exceeded by ω
∗ on submesh (
e
^
) has to be zero, suggesting that J
0
(
e
^
)
=
0. Otherwise, the clamped-end modes can be enlarged infinitely, resulting in the unsatisfaction of the given error tolerance.Take a submesh element(
α
^
i
,
α
^
i
+
1
) on element (e), for instance; its corresponding dynamic stiffness matrix [
k
^
i
]
(
e
^
) 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
(
e
^
) and F
(
e
^
), respectively; the submesh element stiffness matrix k
^
i
(
e
^
) satisfies(17)
k
^
i
e
^
d
e
^
=
F
e
^
.By applying the boundary conditions in (12) in turns again, (17) is expanded as(18)
k
^
i
e
^
B
e
^
=
C
e
^
,where(19)
B
e
^
=
d
j
=
1
e
^
d
j
=
2
e
^
d
j
=
3
e
^
d
j
=
4
e
^
d
j
=
5
e
^
d
j
=
6
e
^
d
j
=
7
e
^
d
j
=
8
e
^
,
C
e
^
=
F
j
=
1
e
^
F
j
=
2
e
^
F
j
=
3
e
^
F
j
=
4
e
^
F
j
=
5
e
^
F
j
=
6
e
^
F
j
=
7
e
^
F
j
=
8
e
^
.As J
0
(
e
^
)
=
0, B
(
e
^
) is nonsingular (otherwise ω
∗ is one of the clamped-end frequencies, which is conflicted with J
0
(
e
^
)
=
0), and (20) is obtained:(20)
k
^
i
e
^
=
C
e
^
B
e
^
-
1
.Since K
(
e
^
)
ω
∗ can be obtained from regular assembly of [
k
^
i
]
(
e
^
), 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 J
0 for nonuniform members. It can be verified that, with the decrease of the length of a shell segment along the meridian, J
0 converges to zero.By employing the W-W algorithm repeatedly, the upper and lower bound of each frequency under a given circumferential wave numbern can be predicted. Recursive Newton’s method is employed to obtain the exact natural frequencies and modes under the specificn, 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.


## 5. 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.J
0, J
K, and J are computed, showing the convergence of the method presented in the previous section.


### 5.1. 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 isd
s
=
R
d
α with α
=
9
0
∘ and R
=
∞. The radius r 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(21)
J
d
z
d
s
=
S
z
.


#### 5.1.1. Convergence Study

Consider a circular cylindrical shell with the lengthL = 24 m, radius r
=
2
m, thickness h
=
0.1
m, Young’s modulus E
=
2.0
×
1
0
11
N
/
m
2, Poisson’s ratio ν
=
0.3, and density ρ
=
7800
k
g
/
m
3. Shell segments are divided equally along the meridian. J
0, J
K, and J with different number of segment membersne at n
=
2, ω
∗
=
3000
r
a
d
/
s are given in Table 1.Table 1
J
0, J
K, and J of a circular cylindrical shell with different ne at n
=
2, ω
∗ = 3000 rad/s.
n
e
J
0
J
K
J
3
24
1
25
4
24
1
25
5
20
5
25
6
18
7
25
8
16
9
25
10
20
5
25
12
12
13
25
16
16
9
25
18
18
7
25
19
0
25
25
20
0
25
25It can be observed that the sum ofJ
0 and J
K always converges to J. With the increase ofne, J
0 will be reduced to zero, suggesting that the calculation of J
0 can be avoided should each shell segment member be short enough. In the following studies,ne is always large enough so that J
0 = 0 can be guaranteed if not specifically stated.


#### 5.1.2. Axisymmetric Vibrations (n
=
0)

For axisymmetric vibration of a Shear Diaphragm-Shear Diaphragm (SD-SD) circular cylindrical shell whose boundary conditions are(22)
v
=
w
=
M
α
=
N
α
=
0
α
=
α
1
,
α
2
,where M
α and N
α 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.Table 2
The lowest twenty natural frequencies (rad/s) of axisymmetric vibration.
Order
Present
Leissa [5]
Order
Present
Leissa [5]
1
0.000167
0.00
11
2413.20
2413.20
2
411.07
411.12
12
2465.41
2465.41
3
660.67
660.67
13
2466.44
2466.70
4
822.15
822.23
14
2490.30
2490.30
5
1233.22
1233.35
15
2505.37
2505.37
6
1304.62
1304.62
16
2516.54
2516.54
7
1644.29
1644.47
17
2526.38
2526.39
8
1885.65
1885.66
18
2536.29
2536.29
9
2055.37
2055.58
19
2547.15
2547.16
10
2268.60
2268.60
20
2559.65
2559.67To investigate the errors between the present method and the results from existing literature, a notationη in (23) is introduced:(23)
η
=
Ω
current
Ω
literature
-
1
,whereΩ
currentandΩ
literatureare 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 d 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 axial-radial coupled frequencies, as all the points in Figure 6 with η
>
10
-
6 belong to the pure-torsion vibrations.Figure 6
Error index of an axisymmetric vibration circular cylindrical shell between the present method and Leissa [5].Denote bym the number of half-waves in the axial direction. It is known that, for eachm, there are two vibration modes: pure-torsion and axial-radial coupled vibration. Figure 7 shows the first three pure-torsion and axial-radial 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 w increases with the increase ofm, displacement u along the axial direction is still dominant for those axial-radial coupled vibrations when m is small.Figure 7
The first three pure-torsion and axial-radial coupled vibration modes atn
=
0.
(a)
Pure-torsion modes
(b)
Axial-radial coupled modes


#### 5.1.3. Vibration of Cantilevered Circular Cylindrical Shells

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: lengthL= 0.5112 m, radius r
=
0.2162
m, and thicknessh= 0.0015 m; the material properties are as follows: Young’s modulusE
=
1.83
×
1
0
11
N
/
m
2, Poisson’s ratioν
=
0.3, and densityρ
=
7492
k
g
/
m
3. The first three natural frequencies with n from 0 to 4 are tabulated and compared with results from Chung [25] in Table 3, reaching very good agreement for most data points.Table 3
The first three natural frequencies (Hz) of a C-F circular cylindrical shell withn from 0 to 4.
n
m
Present
Chung [25]
0
1
1498.95
1497.90
2
2353.11
2351.51
3
3593.53
Not given
1
1
855.92
855.10
2
2320.60
2318.98
3
3076.91
3076.05
2
1
404.70
403.72
2
1438.61
1437.11
3
2489.14
2487.60
3
1
226.04
223.34
2
930.13
928.28
3
1836.96
1834.82
4
1
178.26
171.77
2
647.56
644.48
3
1370.54
1367.64Figure8 gives the first three normalised vibration modes at n
=
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.Figure 8
Normalised modes of the lowest three frequencies of the C-F circular cylindrical shell atn
=
2.


### 5.2. 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 Figure9 with its meridian equation given as(24)
r
2
a
2
-
z
2
b
2
=
1
.Consider a cooling tower shell witha= 25.6032 m,b= 63.9064 m, t
z = 18.5928 m, b
z = 82.1944 m, thicknessh= 0.127 m, density ρ
=
2404.6
k
g
/
m
3, Young’s modulus E
=
2.069
×
1
0
10
N
/
m
2, and Poisson’s ratio ν
=
0.15. Nonuniform shell segment elements are obtained by equally dividing the angle α from α
2 to α
1 in Figure 3(b).Figure 9
A hyperboloid shell of revolution.Table4 presents J
0, J
K, and J of the cooling tower shell with differentne at n
=
0,ω
∗ = 100 rad/s, showing that the present method is well capable of capturing J
0 for nonuniform shell segments.Table 4
J
0, J
K, and J of a cooling tower shell with different ne at n
=
0, ω
∗ = 100 rad/s.
ne
J
0
J
K
J
5
10
4
14
10
8
6
14
15
7
7
14
20
6
8
14
25
4
10
14
30
4
10
14
35
2
12
14
40
2
12
14
45
2
12
14
50
1
13
14
55
0
14
14Table5 lists the first three natural frequencies of the cooling tower shell with n
=
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.Table 5
Natural frequencies (Hz) of a cooling tower shell.
n
m
Present
Deb Nath [26]
0
1
6.2303
6.2329
2
7.7486
7.7521
3
11.417
11.429
1
1
3.2883
3.2897
2
6.7905
6.7933
3
10.518
10.523
2
1
1.7667
1.7664
2
3.6883
3.6893
3
6.9529
6.9556


### 5.3. Spherical Shells

For spherical shells,R is constant and held positive, and r
=
R
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 (r
=
0), a small circular hole ϕ
o
=
0.125
%
ϕ
N is introduced by de Souza and Croll [27]. The same holes are created in the present dynamic stiffness analysis.Figure 10
9
0
∘ and 4
5
∘ spherical shell caps.
(a)
(b)The first four nondimensional natural frequencies (inω
0
=
E
/
ρ
R
2
1
-
ν
2) with n varying from 1 to 3 of the two spherical shell caps are given in Tables 6 and 7, respectively.Table 6
Natural frequencies of a9
0
∘ spherical shell cap (×
ω
0).
m
n
1
2
3
Present
Ref. [27]
Present
Ref. [27]
Present
Ref. [27]
1
0.5374
0.5385
0.8582
0.8580
0.9001
0.9010
2
0.8468
0.8472
0.9172
0.9191
0.9320
0.9324
3
0.9136
0.9141
0.9384
0.9269
0.9488
0.9547
4
0.9372
0.9366
0.9529
0.9577
0.9653
0.9719Table 7
Natural frequencies of a4
5
∘ spherical shell cap (×
ω
0).
m
n
1
2
3
Present
Ref. [27]
Present
Ref. [27]
Present
Ref. [27]
1
0.9208
0.9239
0.9459
0.9462
0.9588
0.9607
2
0.9571
0.9582
0.9858
0.9886
1.0045
1.0098
3
0.9905
0.9975
1.0396
1.0434
1.0695
1.0776
4
1.0548
1.0613
1.1158
1.1232
1.1703
1.1845The 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 Figure11 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 4
5
∘ and 9
0
∘ spherical shell caps and (
2
) errors do not necessarily increase along with the circumferential wave number n or half-wave numberm, suggesting that the accuracy of higher-order frequencies is not affected by the presented dynamic stiffness method.Figure 11
Error index of spherical shell caps between the present method and de Souza and Croll [27].
(a)
9
0
∘shell cap
(b)
4
5
∘shell cap


## 5.1. 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 isd
s
=
R
d
α with α
=
9
0
∘ and R
=
∞. The radius r 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(21)
J
d
z
d
s
=
S
z
.


### 5.1.1. Convergence Study

Consider a circular cylindrical shell with the lengthL = 24 m, radius r
=
2
m, thickness h
=
0.1
m, Young’s modulus E
=
2.0
×
1
0
11
N
/
m
2, Poisson’s ratio ν
=
0.3, and density ρ
=
7800
k
g
/
m
3. Shell segments are divided equally along the meridian. J
0, J
K, and J with different number of segment membersne at n
=
2, ω
∗
=
3000
r
a
d
/
s are given in Table 1.Table 1
J
0, J
K, and J of a circular cylindrical shell with different ne at n
=
2, ω
∗ = 3000 rad/s.
n
e
J
0
J
K
J
3
24
1
25
4
24
1
25
5
20
5
25
6
18
7
25
8
16
9
25
10
20
5
25
12
12
13
25
16
16
9
25
18
18
7
25
19
0
25
25
20
0
25
25It can be observed that the sum ofJ
0 and J
K always converges to J. With the increase ofne, J
0 will be reduced to zero, suggesting that the calculation of J
0 can be avoided should each shell segment member be short enough. In the following studies,ne is always large enough so that J
0 = 0 can be guaranteed if not specifically stated.


### 5.1.2. Axisymmetric Vibrations (n
=
0)

For axisymmetric vibration of a Shear Diaphragm-Shear Diaphragm (SD-SD) circular cylindrical shell whose boundary conditions are(22)
v
=
w
=
M
α
=
N
α
=
0
α
=
α
1
,
α
2
,where M
α and N
α 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.Table 2
The lowest twenty natural frequencies (rad/s) of axisymmetric vibration.
Order
Present
Leissa [5]
Order
Present
Leissa [5]
1
0.000167
0.00
11
2413.20
2413.20
2
411.07
411.12
12
2465.41
2465.41
3
660.67
660.67
13
2466.44
2466.70
4
822.15
822.23
14
2490.30
2490.30
5
1233.22
1233.35
15
2505.37
2505.37
6
1304.62
1304.62
16
2516.54
2516.54
7
1644.29
1644.47
17
2526.38
2526.39
8
1885.65
1885.66
18
2536.29
2536.29
9
2055.37
2055.58
19
2547.15
2547.16
10
2268.60
2268.60
20
2559.65
2559.67To investigate the errors between the present method and the results from existing literature, a notationη in (23) is introduced:(23)
η
=
Ω
current
Ω
literature
-
1
,whereΩ
currentandΩ
literatureare 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 d 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 axial-radial coupled frequencies, as all the points in Figure 6 with η
>
10
-
6 belong to the pure-torsion vibrations.Figure 6
Error index of an axisymmetric vibration circular cylindrical shell between the present method and Leissa [5].Denote bym the number of half-waves in the axial direction. It is known that, for eachm, there are two vibration modes: pure-torsion and axial-radial coupled vibration. Figure 7 shows the first three pure-torsion and axial-radial 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 w increases with the increase ofm, displacement u along the axial direction is still dominant for those axial-radial coupled vibrations when m is small.Figure 7
The first three pure-torsion and axial-radial coupled vibration modes atn
=
0.
(a)
Pure-torsion modes
(b)
Axial-radial coupled modes


### 5.1.3. Vibration of Cantilevered Circular Cylindrical Shells

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: lengthL= 0.5112 m, radius r
=
0.2162
m, and thicknessh= 0.0015 m; the material properties are as follows: Young’s modulusE
=
1.83
×
1
0
11
N
/
m
2, Poisson’s ratioν
=
0.3, and densityρ
=
7492
k
g
/
m
3. The first three natural frequencies with n from 0 to 4 are tabulated and compared with results from Chung [25] in Table 3, reaching very good agreement for most data points.Table 3
The first three natural frequencies (Hz) of a C-F circular cylindrical shell withn from 0 to 4.
n
m
Present
Chung [25]
0
1
1498.95
1497.90
2
2353.11
2351.51
3
3593.53
Not given
1
1
855.92
855.10
2
2320.60
2318.98
3
3076.91
3076.05
2
1
404.70
403.72
2
1438.61
1437.11
3
2489.14
2487.60
3
1
226.04
223.34
2
930.13
928.28
3
1836.96
1834.82
4
1
178.26
171.77
2
647.56
644.48
3
1370.54
1367.64Figure8 gives the first three normalised vibration modes at n
=
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.Figure 8
Normalised modes of the lowest three frequencies of the C-F circular cylindrical shell atn
=
2.


## 5.1.1. Convergence Study

Consider a circular cylindrical shell with the lengthL = 24 m, radius r
=
2
m, thickness h
=
0.1
m, Young’s modulus E
=
2.0
×
1
0
11
N
/
m
2, Poisson’s ratio ν
=
0.3, and density ρ
=
7800
k
g
/
m
3. Shell segments are divided equally along the meridian. J
0, J
K, and J with different number of segment membersne at n
=
2, ω
∗
=
3000
r
a
d
/
s are given in Table 1.Table 1
J
0, J
K, and J of a circular cylindrical shell with different ne at n
=
2, ω
∗ = 3000 rad/s.
n
e
J
0
J
K
J
3
24
1
25
4
24
1
25
5
20
5
25
6
18
7
25
8
16
9
25
10
20
5
25
12
12
13
25
16
16
9
25
18
18
7
25
19
0
25
25
20
0
25
25It can be observed that the sum ofJ
0 and J
K always converges to J. With the increase ofne, J
0 will be reduced to zero, suggesting that the calculation of J
0 can be avoided should each shell segment member be short enough. In the following studies,ne is always large enough so that J
0 = 0 can be guaranteed if not specifically stated.


## 5.1.2. Axisymmetric Vibrations (n
=
0)

For axisymmetric vibration of a Shear Diaphragm-Shear Diaphragm (SD-SD) circular cylindrical shell whose boundary conditions are(22)
v
=
w
=
M
α
=
N
α
=
0
α
=
α
1
,
α
2
,where M
α and N
α 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.Table 2
The lowest twenty natural frequencies (rad/s) of axisymmetric vibration.
Order
Present
Leissa [5]
Order
Present
Leissa [5]
1
0.000167
0.00
11
2413.20
2413.20
2
411.07
411.12
12
2465.41
2465.41
3
660.67
660.67
13
2466.44
2466.70
4
822.15
822.23
14
2490.30
2490.30
5
1233.22
1233.35
15
2505.37
2505.37
6
1304.62
1304.62
16
2516.54
2516.54
7
1644.29
1644.47
17
2526.38
2526.39
8
1885.65
1885.66
18
2536.29
2536.29
9
2055.37
2055.58
19
2547.15
2547.16
10
2268.60
2268.60
20
2559.65
2559.67To investigate the errors between the present method and the results from existing literature, a notationη in (23) is introduced:(23)
η
=
Ω
current
Ω
literature
-
1
,whereΩ
currentandΩ
literatureare 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 d 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 axial-radial coupled frequencies, as all the points in Figure 6 with η
>
10
-
6 belong to the pure-torsion vibrations.Figure 6
Error index of an axisymmetric vibration circular cylindrical shell between the present method and Leissa [5].Denote bym the number of half-waves in the axial direction. It is known that, for eachm, there are two vibration modes: pure-torsion and axial-radial coupled vibration. Figure 7 shows the first three pure-torsion and axial-radial 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 w increases with the increase ofm, displacement u along the axial direction is still dominant for those axial-radial coupled vibrations when m is small.Figure 7
The first three pure-torsion and axial-radial coupled vibration modes atn
=
0.
(a)
Pure-torsion modes
(b)
Axial-radial coupled modes


## 5.1.3. Vibration of Cantilevered Circular Cylindrical Shells

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: lengthL= 0.5112 m, radius r
=
0.2162
m, and thicknessh= 0.0015 m; the material properties are as follows: Young’s modulusE
=
1.83
×
1
0
11
N
/
m
2, Poisson’s ratioν
=
0.3, and densityρ
=
7492
k
g
/
m
3. The first three natural frequencies with n from 0 to 4 are tabulated and compared with results from Chung [25] in Table 3, reaching very good agreement for most data points.Table 3
The first three natural frequencies (Hz) of a C-F circular cylindrical shell withn from 0 to 4.
n
m
Present
Chung [25]
0
1
1498.95
1497.90
2
2353.11
2351.51
3
3593.53
Not given
1
1
855.92
855.10
2
2320.60
2318.98
3
3076.91
3076.05
2
1
404.70
403.72
2
1438.61
1437.11
3
2489.14
2487.60
3
1
226.04
223.34
2
930.13
928.28
3
1836.96
1834.82
4
1
178.26
171.77
2
647.56
644.48
3
1370.54
1367.64Figure8 gives the first three normalised vibration modes at n
=
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.Figure 8
Normalised modes of the lowest three frequencies of the C-F circular cylindrical shell atn
=
2.


## 5.2. 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 Figure9 with its meridian equation given as(24)
r
2
a
2
-
z
2
b
2
=
1
.Consider a cooling tower shell witha= 25.6032 m,b= 63.9064 m, t
z = 18.5928 m, b
z = 82.1944 m, thicknessh= 0.127 m, density ρ
=
2404.6
k
g
/
m
3, Young’s modulus E
=
2.069
×
1
0
10
N
/
m
2, and Poisson’s ratio ν
=
0.15. Nonuniform shell segment elements are obtained by equally dividing the angle α from α
2 to α
1 in Figure 3(b).Figure 9
A hyperboloid shell of revolution.Table4 presents J
0, J
K, and J of the cooling tower shell with differentne at n
=
0,ω
∗ = 100 rad/s, showing that the present method is well capable of capturing J
0 for nonuniform shell segments.Table 4
J
0, J
K, and J of a cooling tower shell with different ne at n
=
0, ω
∗ = 100 rad/s.
ne
J
0
J
K
J
5
10
4
14
10
8
6
14
15
7
7
14
20
6
8
14
25
4
10
14
30
4
10
14
35
2
12
14
40
2
12
14
45
2
12
14
50
1
13
14
55
0
14
14Table5 lists the first three natural frequencies of the cooling tower shell with n
=
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.Table 5
Natural frequencies (Hz) of a cooling tower shell.
n
m
Present
Deb Nath [26]
0
1
6.2303
6.2329
2
7.7486
7.7521
3
11.417
11.429
1
1
3.2883
3.2897
2
6.7905
6.7933
3
10.518
10.523
2
1
1.7667
1.7664
2
3.6883
3.6893
3
6.9529
6.9556


## 5.3. Spherical Shells

For spherical shells,R is constant and held positive, and r
=
R
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 (r
=
0), a small circular hole ϕ
o
=
0.125
%
ϕ
N is introduced by de Souza and Croll [27]. The same holes are created in the present dynamic stiffness analysis.Figure 10
9
0
∘ and 4
5
∘ spherical shell caps.
(a)
(b)The first four nondimensional natural frequencies (inω
0
=
E
/
ρ
R
2
1
-
ν
2) with n varying from 1 to 3 of the two spherical shell caps are given in Tables 6 and 7, respectively.Table 6
Natural frequencies of a9
0
∘ spherical shell cap (×
ω
0).
m
n
1
2
3
Present
Ref. [27]
Present
Ref. [27]
Present
Ref. [27]
1
0.5374
0.5385
0.8582
0.8580
0.9001
0.9010
2
0.8468
0.8472
0.9172
0.9191
0.9320
0.9324
3
0.9136
0.9141
0.9384
0.9269
0.9488
0.9547
4
0.9372
0.9366
0.9529
0.9577
0.9653
0.9719Table 7
Natural frequencies of a4
5
∘ spherical shell cap (×
ω
0).
m
n
1
2
3
Present
Ref. [27]
Present
Ref. [27]
Present
Ref. [27]
1
0.9208
0.9239
0.9459
0.9462
0.9588
0.9607
2
0.9571
0.9582
0.9858
0.9886
1.0045
1.0098
3
0.9905
0.9975
1.0396
1.0434
1.0695
1.0776
4
1.0548
1.0613
1.1158
1.1232
1.1703
1.1845The 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 Figure11 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 4
5
∘ and 9
0
∘ spherical shell caps and (
2
) errors do not necessarily increase along with the circumferential wave number n or half-wave numberm, suggesting that the accuracy of higher-order frequencies is not affected by the presented dynamic stiffness method.Figure 11
Error index of spherical shell caps between the present method and de Souza and Croll [27].
(a)
9
0
∘shell cap
(b)
4
5
∘shell cap


## 6. 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 one-dimensional problems with respect to different circumferential wave numbern, 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 J
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.



---

*Source: 4513520-2016-09-07.xml*