Open access

Vibrations of Cylindrical Shells

Written By

Tiejun Yang, Wen L. Li and Lu Dai

Submitted: 07 March 2012 Published: 02 October 2012

DOI: 10.5772/51816

Chapter metrics overview

2,914 Chapter Downloads

View Full Metrics

1. Introduction

Beams, plates and shells are the most commonly-used structural components in industrial applications. In comparison with beams and plates, shells usually exhibit more complicated dynamic behaviours because the curvature will effectively couple the flexural and in-plane deformations together as manifested in the fact that all three displacement components simultaneously appear in each of the governing differential equations and boundary conditions. Thus it is understandable that the axial constraints can have direct effects on a predominantly radial mode. For instance, it has been shown that the natural frequencies for the circumferential modes of a simply supported shell can be noticeably modified by the constraints applied in the axial direction [1]. Vibrations of shells have been extensively studied for several decades, resulting in numerous shell theories or formulations to account for the various effects associated with deformations or stress components.

Expressions for the natural frequencies and modes shapes can be derived for the classical homogeneous boundary conditions [2-9]. Wave propagation approach was employed by several researchers [10-13] to predict the natural frequencies for finite circular cylindrical shells with different boundary conditions. Because of the complexity and tediousness of the (exact) solution procedures, approximate procedures such as the Rayleigh-Ritz methods or equivalent energy methods have been widely used for solving shell problems [14-18]. In the Rayleigh-Ritz methods, the characteristic functions for a “similar” beam problem are typically used to represent all three displacement components, leading to a characteristic equation in the form of cubic polynomials. Assuming that the circumferential wave length is smaller than the axial wave length, Yu [6] derived a simple formula for calculating the natural frequencies directly from the shell parameters and the frequency parameters for the analogous beam case. Soedel [19] improved and generalized Yu’s result by eliminating the short circumferential wave length restriction. However, since the wavenumbers for axial modal function are obtained from beam functions which do not exactly satisfy shell boundary conditions, it is mathematically difficult to access or ensure the accuracy and convergence of such a solution.

The free vibration of shells with elastic supports was studied by Loveday and Rogers [20] using a general analysis procedure originally presented by Warburton [3]. The effect of flexibility in boundary conditions on the natural frequencies of two (lower order) circumferential modes was investigated for a range of restraining stiffness values. The vibrations of circular cylindrical shells with non-uniform boundary constraints were studied by Amabili and Garziera [21] using the artificial spring method in which the modes for the corresponding less-restrained problem were used to expand the displacement solutions. The non-uniform spring stiffness distributions were systematically represented by cosine series and their presence was accounted for in terms of maximum potential energies stored in the springs.

A large number of studies are available in the literature for the vibrations of shells under different boundary conditions or with various complicating features. A comprehensive review of early investigations can be found in Leissa’s book [22]. Some recent progresses have been reviewed by Qatu [23]. Regardless of whether an approximate or an exact solution procedure is employed, the corresponding formulations and implementations usually have to be modified or customized for different boundary conditions. This shall not be considered a trivial task in view that there exist 136 different combinations even considering the simplest (homogeneous) boundary conditions. Thus, it is useful to develop a solution method that can be generally applied to a wide range of boundary conditions with no need of modifying solution algorithms and procedures. Mathematically, elastic supports represent a general form of boundary conditions from which all the classical boundary conditions can be readily derived by simply setting each of the spring stiffnesses to either zero or infinity. This chapter will be devoted to developing a general analytical method for solving shell problems involving general elastically restrained ends.


2. Basic equations and solution procedures

Figure 1 shows an elastically restrained circular cylindrical shell of radius R, thickness h and length L. Each of the eight sets of elastic restraints shall be understood as a distributed spring along the circumference. Let u, v, and w denote the displacements in the axial x, circumferential θ and radial r directions, respectively. The equations of the motions for the shell can be written as

N1x+N12Rθ=ρh2ut2 (a)N12x+N2Rθ=ρh2vt2 (b)Q1x+Q1RθN2R=ρh2wt2 (c)E1

where ρ is the mass density of the shell material, and N1, N12, N, Q1 and Q denote the resultant forces acting on the mid-surface.

Figure 1.

A circular cylindrical shell elastically restrained along all edges.

In terms of the shell displacements, the force and moment components can be expressed as

N1=K(ux+σvRθ+σwR) (a)N2=K(vRθ+wR+σux) (b)N12=K(1σ)2(uRθ+vx) (c)M1=D(2wx2+σ2wR2θ2) (d)M2=D(2wR2θ2+σ2wx2) (e)M12=D(1σ)2wRxθ (f)Q1=M1x+M12Rθ=D(3wx3+3wR2xθ2) (g)Q2=M2Rθ+M12x=D(3wR3θ3+3wRx2θ) (h)E2


K=Eh/(1σ2), (a)D=Eh3/12(1σ2) (b)κ=D/K=h2/12 (c)E3

E and σ are respectively the Young’s modulus and Poisson ratio of the material; M1, M2 and M12 are the bending and twisting moments.

The boundary conditions for an elastically restrained shell can specified as:

at x=0,

Nxk1u=0 (a)Nxθk2v=0 (b)Qx+MxθRθk3w=0 (c)Mx+k4wx=0 (d)E4

at x=L,

Nx+k5u=0 (a)       Nxθ+k6v=0 (b)Qx+MxθRθ+k7w=0 (c)Mxk8wx=0 (d)E5

where k1, k2, …, k8 are the stiffnesses for the restraining springs. The elastic supports represent a set of general boundary conditions, and all the classical boundary conditions can be considered as the special cases when the stiffness for each spring is equal to either zero or infinity.

The above equations are usually referred to as Donnell-Mushtari equations. Flügge’s theory is also widely used to describe vibrations of shells. In terms of the shell displacements, the corresponding force and moment components are written as

N1=K(ux+σvRθ)+KσwRDR2wx2 (a)N12=1σ2K(uRθ+vx)+DR(1σ)2(vRx2wRθx) (b)N2=K(vRθ+wR+σux)+DR(2wR2θ2+wR2) (c)N21=1σ2K(uRθ+vx)+DR(1σ)2(uR2θ+2wRθx) (d)M1=D(2wx2+σ2wR2θ2uRxσvR2θ) (e)M12=D(1σ)(2wRθxvRx) (f)M2=D(wR2+2wR2θ2+σ2wx2) (g)M21=D(1σ)(2wRθx+12uR2θ12vRx) (h)Q1=D(3wx3+1R23wθ2x2uRx2+1σ2R32uθ21+σ2R2vxθ) (i)E6

A shell problem can be solved either exactly or approximately. An exact solution usually implies that both the governing equations and the boundary conditions are simultaneously satisfied exactly on a point-wise basis. Otherwise, a solution is considered approximate in which one or more of the governing equations and boundary conditions are enforced only in an approximate sense. Both solution strategies will be used below.

2.1. An approximate solution based on the Rayleigh-Ritz procedure

Approximate methods based on energy methods or the Rayleigh-Ritz procedures are widely used for the vibration analysis of shells with various boundary conditions and/or complicating factors. In such an approach, the displacement functions are usually expressed as

u(x,θ)=m=0amφum(x) cosnθ (a)v(x,θ)=m=0bmφvm(x) sinnθ (b)w(x,θ)=m=0cmφwm(x) cosnθ (c)E7

where φαm(x), α=u,v, and w, are the characteristic functions for beams with similar boundary conditions. Although characteristic functions generally exist in the forms of trigonometric and hyperbolic functions, they also include some integration and frequency constants that have to be determined from boundary conditions. Consequently, each boundary condition basically calls for a special set of modal data. In the literature the modal parameters are well established for the simplest homogeneous boundary conditions. However, the determination of modal properties for the more complicated elastic boundary supports can become, at least, a tedious task since they have to be re-calculated each time when any of the stiffness values is changed. It should also be noted that calculating the modal properties will typically involve seeking the roots of a nonlinear transcendental equation, which mathematically requires an iterative root searching scheme and careful numerical implementations to ensure no missing data. To overcome this problem, a unified representation of the shell solutions will be adopted here in which the displacements, regardless of boundary conditions, will be invariably expressed as

u(x,θ)=(m=0amcosλmx +pu(x)) cosnθ, (λm=mπL) (a)v(x,θ)=(m=0bmcosλmx +pv(x)) sinnθ ,  (b)w(x,θ)=(m=0cmcosλmx +pw(x)) cosnθ (c)


where pα(x), α=u,v, and w, denote three auxiliary polynomials which satisfy

pu(x)x|x=0=u(x,0)x|x=0=β1, (a)pu(x)x|x=L=u(x,0)x|x=L=β2 (b)pv(x)x|x=0=v(x,π/2)x|x=0=β3, (c)pv(x)x|x=L=v(x,π/2)x|x=L=β4, (d)pw(x)x|x=0=w(x,0)x|x=0=β5, (e)pw(x)x|x=L=w(x,0)x|x=L=β6, (f)3pw(x)x3|x=0=3w(x,0)x3|x=0=β7, (g)3pw(x)x3|x=L=3w(x,0)x3|x=L=β8 (h) E9

It is clear from Eqs. (9) that these auxiliary polynomials are only dependent on the first and third derivatives βi, (i=1,2,…,8) of the displacement solutions on the boundaries. In terms of boundary derivatives, the lowest-order polynomials can be explicitly expressed as [24, 25]

pu(x)=ζ1(x)β1+ζ2(x)β2 (a)pv(x)=ζ1(x)β3+ζ2(x)β4 (b)pw(x)=ζ1(x)β5+ζ2(x)β6+ζ3(x)β7+ζ4(x)β8 (c)E10


ζ(x)T={ζ1(x)ζ2(x) ζ3(x)ζ4(x)}={(6Lx2L23x2)/6L(3x2L2)/6L (15x460Lx3+60L2x28L4)/360L (15x430L2x2+7L4)/360L}E11

This alternative form of Fourier series recognizes the fact that the conventional Fourier series for a sufficiently smooth function f(x) defined on a compact interval [0, L] generally fails to converge at the end points. Introducing the auxiliary functions will ensure the cosine series in Eqs. (8) to converge uniformly and polynomially over the interval, including the end points. As a matter of fact, the polynomial subtraction techniques have been employed by mathematicians as a means to accelerate the convergence of the Fourier series expansion for an explicitly defined function [26-28].

The coefficients βi represent the values of the first and third derivatives of the displacements at the boundaries, and are hence related to the unknown Fourier coefficients for the trigonometric terms. The relationships between the constants and the expansion coefficients can be derived either exactly or approximately.

In seeking an approximate solution based on an energy method, the solution is not required to explicitly satisfy the force or natural boundary conditions. Accordingly, the derivative parameters βi in Eqs. (10) will be here determined from a simplified set of the boundary conditions, that is,

uxk^1,5u=0 (a)(1μ)2vxk^2,6v=0 (b)κ3wx3±k^3,7w=0 (c)κ2wx2k^4,8wx=0 (d)E12



By substituting Eqs. (8) and (10) into (12), one will obtain

{β1β2}=m=0Hu1Qumam (a){β3β4}=m=0Hv1Qvmbm (b){ β5β6β7β8 }T=m=0Hw1Qwmcm (c)E14


Hu=[k^1L3+1k^1L6k^5L6k^5L3+1] (a)Qum={ k^1(1)m+1k^5}T (b)Hv=[k^2L3+1μ2k^2L6k^6L6k^6L3+1μ2] (c)Qvm={ k^2(1)m+1k^6}T (d)Hw=[k^3L3k^3L6k^3L345+κ7k^3L3360k^7L6k^7L37k^7L3360k^7L345+κk^4+κLκLκL3κL6κLk^8+κLκL6κL3] (e)E15



In light of Eqs. (13), Eqs (8) can be reduced to Eqs. (7) with the axial functions being defined as



Q¯wm= Q˜wm   (a)Q¯αm= {Q˜αm00}      (b)Q˜αm=(Hα)-1Qαm  (c)E18

Since the boundary conditions are not exactly satisfied by the displacements such constructed, the Rayleigh-Ritz procedure will be employed to find a weak form of solution. The current solution is noticeably different from the conventional Rayleigh-Ritz solutions in that: a) the shell displacements are expressed in terms of three independent sets of axial functions, rather than a single (set of) beam function(s), b) the basis functions in each displacement expansion constitutes a complete set so that the convergence of the Rayleigh-Ritz solution is guaranteed mathematically, and c) it does not suffer from the well-known numerical instability problem related to the higher order beam functions or polynomials. More importantly, the current method is that it provides a unified solution to a wide variety of boundary conditions.

The potential energy consistent with the Donnell-Mushtari theory can be expressed from

V=K20L02π{(ux+vRθ+wR)22(1ν)ux(vRθ+wR)+(1ν)2(vx+uRθ)2+κ[(2wx2+2wR2θ2)22(1ν)(2wx22wR2θ2(2wRxθ)2)]}Rdxdθ+1/202π[k1u2+k2v2+k3w2+k4(w/x)2]x=0Rdθ+1/202π[k5u2+k6v2+k7w2+k8(w/x)2]x=LRdθ (a)and the total kinetic energy is calculated fromT=120L02πρh[(u/t)2+(v/t)2+(w/t)2] Rdxdθ (b)E19

By minimizing the Lagrangian L=V-T against all the unknown expansion coefficients, a final system of linear algebraic equations can be derived



a¯={ a00,a01,...,amn,..}T,    (a)b¯={ b01,b02,...,bmn,...}T,    (b)c¯={ c00,c01,...,cmn,...}T,    (c)E21
Λmn,m'n'ss=δnn'[Iuu,11mm'+(1μ)n22R2Iuu,00mm'+2Lk^1φum(0)φum'(0)+2Lk^5φum(L)φum'(L)] (a)Λmn,m'n'sθ=δnn'[μnRIuv,10mm'(1μ)n2RIuv,01mm'] (b)Λmn,m'n'sr=δnn'μRIuw,10mm' (c)Λmn,m'n'θθ=δnn'[n2R2Ivv,00mm'+(1μ)2Ivv,11mm'+2Lk^2φvm(0)φvm'(0)+2Lk^6φvm(L)φvm'(L)] (d)Λmn,m'n'θr=δnn'nR2Ivw,00mm' (e)Λmn,m'n'θθ=δnn'{1R2Iww,00mm'+κ[Iww,22mm'+n4R4Iww,00mm'+2(1μ)n2R2Iww,11mm'                     μn2R2(Iww,02mm'+Iww,20mm')]+2Lk^3φwm(0)φwm'(0)+2Lk^7φwm(L)φwm'(L)+ (f)                     +2Lk^4φwm(0)xφwm'(0)x+2Lk^8φwm(L)xφwm'(L)x}E22
Mmn,m'n'ss=δnn'ρhIuu,00mm' (a)Mmn,m'n'θθ=δnn'ρhIvv,00mm' (b)Mmn,m'n'rr=δnn'ρhIww,00mm' (c)E23


Iαβ,pqmm'=2/L0Lpϕαmxpqϕβm'xqdx(α,β=u,v,w) E24

The integrals in Eq. (23) can be calculated analytically; for instance

Iαβ,00mm'=2 L0Lφαmφβm' dx 
=2 L0L(cosλmx+ζ(x)TQ¯αm)(cosλm'x+ζ(x)TQ¯βm') dx=2L0L(cosλmxcosλm'x+cosλm'xζ(x)TQ¯αm+cosλmxζ(x)TQ¯βm'+Q¯αmTζ(x)ζ(x)TQ¯βm') dx =εmδmm'+Sαm'm+Sβmm'+Zαβmm'E25


 εm=(1+δm0) (a)Sαm'm=Pcm'Q¯αm (b)Ζαβmm'=Q¯αmTΞQ¯βm'  (c)E26


Pcm=2L0Lζw(x)Tcosλmx dxE27
=2L{{ 0000}T,                                           for m=0{1λm2(1)mλm21λm4(1)m+1λm4}T,              for m0 (a)Ξ=2/L0Lζ(x)Tζ(x)dx=[2L2457L21802L245sym.4L494531L475602L6472531L475604L4945127L63024002L64725] (b)  E28

2.2. A strong form of solution based on Flügge’s equations

As aforementioned, the displacement expressions in terms of beam functions cannot exactly satisfy the shell boundary conditions; instead they are made to satisfy the boundary conditions in a weak sense via the use of the Rayleigh-Ritz procedure. To overcome this problem, the displacement expressions, Eqs. (8), will now be generalized to

u(x,θ)=n=0(m=0Amncosλmx+pnu(x))cosnθ (a)v(x,θ)=n=0(m=0Bmncosλmx+pnv(x))sinnθ (b)w(x,θ)=n=0(m=0Cmncosλmx+pnw(x))cosnθ (c)E29

which represent a 2-D version of the improved Fourier series expansions, Eqs. (8).

To demonstrate the flexibility in choosing the auxiliary functionspnu(x), pnv(x)andpnw(x), an alternative set is used below:

pnu(x)=Λnuα(x) (a)pnv(x)=Λnvα(x) (b)pnw(x)=Λnwβ(x) (c)E30

hereΛnu=[an bn], Λnv=[cn dn], Λnw=[en fn gn hn]with an, bn,..., gn and hn being the unknown coefficients to be determined; α(x)={α1(x)α2(x)}Tand β(x)={β1(x)β2(x)β3(x)β4(x)}T and with their elements being defined as

α1(x)=x(xl1)2 (a)α2(x)=xl2(xl1) (b)E31


β1(x)=9l4πsin(πx2l)l12πsin(3πx2l) (a)β2(x)=9l4πcos(πx2l)l12πcos(3πx2l) (b)β3(x)=l3π3sin(πx2l)l33π3sin(3πx2l) (c)β4(x)=l3π3cos(πx2l)l33π3cos(3πx2l) (d)E32

In Eqs. (27), the sums of x-related terms are here understood as the series expansions in x-direction, rather than characteristic functions for a beam with “similar” boundary condition. This distinction is important in that the boundary conditions and governing differential equations can now be exactly satisfied on a point-wise basis; that is, the solution can be found in strong form, as described below.

Substituting Eqs. (6) and (27) into (4) and (5) will lead to

an(7σl3πR+3πRγ4l)fn(4σl33π3R+lγRπ)hn =k1Km=0AmnσnRm=0Bmnm=0(σR+λm2γR)Cmn (a)1σ2(1+γ)cn+(1σ)γn2en=(1σ)n2Rm=0Amn+k2Km=0Bmn (b)4lRan2lRbn+(3σ)n2R2cn+(2σ)n2R2en+7lk33πDfngn+4l3k33π3Dhn=m=0(λm2R(1σ)n22R3)Amn+k3Dm=0Cmn (c)1Ran+(3π4l+7νln23πR2)fn+(lπ+4νl3n23π3R2)hnk4Den=σnR2m=0Bmn+m=0(λm2+σn2R2)Cmn (d)bn(3πγR4l+7σl3πR)en(lγRπ+4σl33π3R)gn=k5Km=0cos(mπ)Amn+nσRm=0cos(mπ)Bmn+m=0(σR+λm2γR)cos(mπ)Cmn (e)1σ2(1+γ)dn(1σ)nγ2fn=(1σ)n2Rm=0cos(mπ)Amn+k6Km=0cos(mπ)Bmn (f)2lRan4lRbn(3σ)n2R2dn7lk73πDen(2σ)n2R2fn4l3k73π3Dgn+hn=m=0(cos(mπ)λm2R(1σ)cos(mπ)n22R3)Amn+k7Dm=0cos(mπ)Cmn (g)1Rbn+(3π4l+7σln23πR2)enk8Dfn+(lπ+4σl3n23π3R2)gn=σnR2m=0cos(mπ)Bmnm=0(λm2+σn2R2)cos(mπ)Cmn (h)E33

Equations (31) represent a set of constraint conditions between the unknown (boundary) constants, an, bn,..., gn and hn, and the Fourier expansion coefficients Amn, Bmn, and Cmn ( m, n = 0, 1, 2,... ). The constraint equations (31a-h) can be rewritten more concisely, in matrix form, as


The elements of the coefficient matrices can be readily derived from Eqs. (31); for example, Eq. (31a) implies

{L31}n,n=δnn (a){L36}n,n=(7σl3πR+3πγR4l)δnn (b){L38}n,n=(4σl33π3R+lγRπ)δnn (c){L32}n,n={L33}n,n={L34}n,n={L35}n,n={L37}n,n=0(d){S11}mn,n=k1δnn/K (e){S12}mn,n=σnδnn/R (f){S13}mn,n=δnn(σR+m2π2γRl2) (g)E35

Other sub-matrices can be similarly obtained from the remaining equations in Eqs. (31).

In actual numerical calculations, all the series expansions will have to be truncated to m=M and n=N. Thus there is a total number of (M+1)(3N+2)+8N+6 unknown expansion coefficients in the displacement functions. Since Eq. (33) represents a set of 8N+6 equations, additional (M+1)(3N+2) equations are needed to be able to solve for all the unknown coefficients. Accordingly, we will turn to the governing differential equations.

In Flügge’s theory, the equations of motion are given as

N1x+N21Rθ=ρh2ut2 (a)N12x+N2Rθ+(M2R2θ+M12Rx)=ρh2vt2 (b)2M1x2+2M12Rxθ+2M21Rxθ+2M2R2θ2N2R=ρh2wt2 (c)E36

Substituting Eqs. (6) and (37) into Eqs. (34) results in

m=0n=0(λm2(1σ)(1+γ)n22R2)cosλmxAmn+n=0Λnu(α(x)(1σ)(1+γ)n22R2)+(1+σ)2R(m=0n=0λmnsinλmxBmn+n=0nΛnvα(x))m=0n=0(σRλm+γRλm3(1σ)γn22Rλm)sinλmxCmn+n=0ΛnwΤ(σRβ(x)γRβ(x)(1σ)γn22Rβ(x))+ω2ρhK(m=0n=0cosλmxAmn+n=0Λnuα(x))=0 (a)+(1+σ)2R(m=0n=0λmnsinλmxAmnn=0nΛnuα(x))m=0n=0(n2R2+(1σ)(1+3γ)2λm2)cosλmxBmnn=0ΛnvΤ(n2R2α(x)(1σ)(1+3γ)2α(x))m=0n=0(nR2+(3σ)γn2λm2)cosλmxCmnn=0ΛnwΤ(nR2β(x)(3σ)γn2β(x))+ω2ρhK(m=0n=0cosλmxBmn+n=0Λnvα(x))=0 (b)m=0n=0(σR(1σ)γn22R+γRλm2)λmsinλmxAmn+n=0ΛnuΤ(σRα(x)(1σ)γn22Rα(x)γRα(x))+m=0n=0(nR2+(3σ)γn2λm2)cosλmxBmn+n=0ΛnuΤ(nR2α(x)(3σ)γn2α(x))+m=0n=0(1+γ(n2+1)2R2+R2γλm4+2γn2λm2)cosλmxCmn+n=0ΛnwΤ(1+γ(n2+1)2R2β(x)2γn2β(x)+γR2β(x))ω2ρhK(m=0n=0cosλmxCmn+n=0Λnwβ(x))=0 (c)E37

By expanding all non-cosine terms into Fourier cosine series and comparing the like terms, the following matrix equation can be obtained


where E, F, P and Q are coefficient matrices whose elements are given as:

{E11}mn,mn=(λm2+(1σ)(1+γ)n22R2)δmmδnn (a){E12}mn,mn=(1+σ)mnπ2lRχmmδmmδnn (b){E13}mn,mn=(σmπlR+γ(Rλm3(1σ)mn2π2lR))χmmδmmδnn (c){E21}mn,mn=(1+σ)mnπ2lRχmmδmmδnn (d){E22}mn,mn=(n2R2+(1ν)(1+3γ)λm22)δmmδnn (e){E23}mn,mn=(nR2+(3ν)γnλm22)δmmδnn (f){E31}mn,mn=(νλmR+γ((1σ)λmn22RRλm3))χmmδmmδnn (g){E32}mn,mn=(nR2+(3σ)λm2nγ2)δmmδnn (h){E33}mn,mn=(1R2+γ((Rλm+n21R)2+2λm))δmmδnn (i)E39
{F11a}mn,n=(ψ1m(1σ)(1+γ)n22R2ϕ1m)δnn (a){F11b}mn,n=(ψ2m(1σ)(1+γ)n22R2ϕ1m)δnn (b){F12c}mn,n=(1+σ)n2Rφ1mδnn (c){F12d}mn,n=(1+σ)n2Rφ2mδnn (d){F13e}mn,n=((9σ8R+γ(9Rπ232l29(1σ)n216R))κ2m(σ8R+γ(9Rπ232l2(1σ)n216R))κ4m)δnn (e){F13f}mn,n=((9σ8R+γ(9Rπ232l29(1σ)n216R))κ1m+(σ8R+γ(9Rπ232l2(1σ)n216R))κ3m)δnn (f){F13g}mn,n=((σl22π2R+γ(R8(1σ)l2n24π2R))κ2m(σl22π2R+γ(9R8(1σ)l2n24π2R))κ4m)δnn (g){F13h}mn,n=((σl22π2R+γ(R8(1σ)l2n24π2R))κ1m+(σl22π2R+γ(9R8(1σ)l2n24π2R))κ3m)δnn (h){F21a}mn,n=(1+σ)n2Rφ1mδnn (i){F21b}mn,n=(1+σ)n2Rφ2mδnn (j){F22c}mn,n=(n2R2ϕ1m(1σ)(1+3γ)2ψ1m)δnn (k){F22d}mn,n=(n2R2ϕ2m(1σ)(1+3γ)2ψ2m)δnn (l){F23e}mn,n=((9nl4πR2+9π(3σ)nγ32l)κ1m(nl12πR2+3π(3σ)nγ32l)κ3m)δnn (m){F23f}mn,n=((9nl4πR2+9π(3σ)nγ32l)κ2m+(nl12πR2+3π(3σ)nγ32l)κ4m)δnn (n){F23g}mn,n=((nl3π3R2+(3σ)nl8π)κ1m(nl33π3R2+3(3σ)nlγ8π)κ3m)δnn (o){F23h}mn,n=((nl3π3R2+(3σ)nl8π)κ2m+(nl33π3R2+3(3σ)nlγ8π)κ4m)δnn (p){F31a}mn,n=((σR(1σ)γn22R)φ1m6Rγl2δm0)δnn (o){F31b}mn,n=((σR(1σ)γn22R)φ2m6Rγl2δm0)δnn (r){F32c}mn,n=(nR2ϕ1m(3σ)nγ2ψ1m)δnn (s){F32d}mn,n=(nR2ϕ21m(3σ)nγ2ψ2m)δnn (t){F33e}mn,n=(9l4π(1R2+γ((Rπ24l2+n21R)2+π22l2))κ1ml3π(14R2+γ((9Rπ28l2+n212R)2+9π28l2))κ3m)δnn (u){F33f}mn,n=(9l4π(1R2+γ((Rπ24l2+n21R)2+π22l2))κ2ml3π(14R2+γ((9Rπ28l2+n212R)2+9π28l2))κ4m)δnn (v){F33g}mn,n=(lπ(l2π2R2+γ((Rπ4l+(n21)lπR)212))κ1ml3π(l2π2R2+γ((9Rπ4l+(n21)lπR)292))κ3m)δnn (w){F33h}mn,n=(lπ(l2π2R2+γ((Rπ4l+(n21)lπR)212))κ2ml3π(l2π2R2+γDesign Science, Inc.((9Rπ4l+(n21)lπR)292))κ4m)δnn (x)E40
{P11}mn,mn={P22}mn,mn={P33}mn,mn=δmmδnn (a){Q11a}mn,n={Q22c}mn,n=ϕ1mδnn (b){Q11b}mn,n={Q22d}mn,n=ϕ2mδnn (c){Q33e}mn,n=l4π(9κ1mκ3m3)δnn (d){Q33f}mn,n=l4π(9κ2m+κ4m3)δnn (e){Q33g}mn,n=l3π3(κ1mκ3m3)δnn (f){Q33h}mn,n=l3π3(κ2m+κ4m3)δnn (g)E41

The symbolsκ1m, κ2m, κ3m, κ4m, ϕ1m, ϕ2m, φ1m, φ2m, ψ1m, ψ2m, and χmi in the above equations are defined as

sin(π2lx)=m=0κ1mcosλmx; κ1m={       2π             m=0,4(14m2)π      m0, (a,b) cos(π2lx)=m=0κ2mcosλmx; κ2m={       2π              m=0,4(1)m(14m2)π       m0, (c,d)sin(3π2lx)=m=0κ3mcosλmx; κ3m={      23π            m=0,12(94m2)π      m0, (e,f)cos(3π2lx)=m=0κ4mcosλmx; κ4m={   23π             m=0,12(1)m+1(94m2)π      m0, (g,h)α1(x)=m=0ϕ1mcosλmx; ϕ1m={                  l12                   m=0,2l(m2π26+6(1)m)m4π4   m0, (i,j)α2(x)=m=0ϕ2mcosλmx; ϕ2m={                 l12                      m=0,2l(m2π2(1)m+66(1)m)m4π4   m0, (k,l)α1(x)=m=0φ1mcosλmx; φ1m={          0             m=0,4(2+(1)m)m2π2    m0, (m,n)α2(x)=m=0φ2mcosλmx; φ2m={           0              m=0,4(1+2(1)m)m2π2    m0, (o,p)α1(x)=m=0ψ1mcosλmx; ψ1m={          1l              m=0,12(1+(1)m)lm2π2     m0, (q,r)α2(x)=m=0ψ2mcosλmx; ψ2m={           1l               m=0,12(1+(1)m)lm2π2    m0, (s,t)sinλmx=iχimcosλix=sinλix=mχmicosλmx; χmi={                    0                         i=0,{    1(1)iiπ          m=0,2i[(1)m+i1](m2i2)π     m0,      i0. (u,v)E42

All the unmentioned elements in matrices P and Q are identically equal to zero.

Equations (32) and (36) can be combined into


where K=E+FL-1S andM=P+QL-1S.

The final system of equations, Eq. (19) or (41), represents a standard characteristic equation for a matrix eigen-problem from which all the eigenvalues and eigenvectors can be readily calculated. It should be mentioned that the elements in each eigenvector are actually the expansion coefficients for the corresponding mode; its “physical” mode shape can be directly obtained from Eqs. (7) or (27).

In the above discussions, the stiffness distribution for each restraining spring is assumed to be axisymmetric or uniform along the circumference. However, this restriction is not necessary. For non-uniform elastic boundary restraints, the displacement expansions, Eq. (27), shall be used, and any and all of stiffness constants can be simply understood as varying with spatial angle θ. For simplicity, we can universally expand these functions into standard cosine series and modify Eq. (31) accordingly to reflect this complicating factor.


3. Results and discussion

Several numerical examples will be given below to verify the two solution strategies described earlier.

3.1. Results about the approximate Rayleigh-Ritz solution

We first consider a familiar simply-supported cylindrical shell. The simply supported boundary condition, Nx=Mx=v=w=0at each end, can be considered as a special case when k2,6=k3,7= and k1,5=k4,8=0 (in actual calculations, infinity is represented by a sufficiently large number). To examine the convergence of the solution, Table 1 shows the frequency parameters, Ω=ωRρ(1σ2)/E, calculated using different numbers of terms in the series expansions. It is seen that the solution converges nicely with only a small number of terms. In the following calculations, the expansions in axial direction will be simply truncated to M=15. Given in Table 2 are the frequencies parameters for some lower-order modes. Exact solution is available for the simply supported case and the results are also shown there for comparison. An excellent agreement is observed between these two sets of results. Although the simply supported boundary condition represents the simplest case in shell analysis, this problem is not trivial in testing the reliability and sophistication of the current solution method. From numerical analysis standpoint, it may actually represent a quite challenging case because of the extreme stiffness values involved. The non-trivialness can also been seen mathematically from the fact that the simple sine function (in the axial direction) in the exact solution is actually expanded as a cosine series expansion in the current solution.

Number of terms used in the series Ω=ωRρ(1σ2)/E
n=0 n=1 n=2 n=3 n=4

Table 1.

Frequency parameters, Ω=ωRρ(1σ2)/E, obtained using different numbers of terms in the displacement expansions.

Mode Ω=ωRρ(1σ2)/E
n=0 n=1 n=2 n=3 n=4
m=1, Current
m=2, Current
m=3, Current

Table 2.

Frequency parameters, Ω=ωRρ(1σ2)/E, for a simply-supported shell; L=4R, h/R=0.05 and μ=0.3.

Next, consider a cylindrical shell clamped at each end, that is,u=v=w=w/x=0. The clamped-clamped boundary condition is a case when the stiffnesses of the restraining springs all become infinitely large. The related shell and material parameters are as follows: L=0.502 m, R=0.0635 m, h=0.00163 m, E=2.1×1011, μ=0.28, and ρ=7800. Listed in Table 3 are some of the lowest natural frequencies for this clamped-clamped shell. The reference results given there are calculated from


whereΩ=ωRρ(1σ2)/E, and the coefficients A0, A1 and A2 are the functions of the modal indices, shell parameters, and the boundary conditions [27]. Equation (42) can be derived from the Rayleigh-Ritz procedure by adopting the beam characteristic functions as the axial functions for all three displacement components. A noticeable difference between these two sets of results may be attributed to the fact that: a) Eq. (42) given in ref. [29] is based on the Flügge shell theory, rather than the Donnell-Mushtari theory, and b) Eq. (42) uses only a single beam characteristic function in contrast to the three complete sets (of basis functions) in the current method.

Mode Current
Eq. (42) Current
Eq. (42)

Table 3.

The natural frequencies in Hz for a clamped-clamped shell; L=0.502 m, R=0.0635 m, h=0.00163 m, E=2.1E+11, μ=0.28, ρ=7800 kg/m3.

Another classical example involves a completely free shell. Vibration of a free-free shell is of particular interest as manifested in the debate between two legendary figures, Rayleigh and Love, about the validity of the inextensional theory of shells. The lower-order modes are typically related to the rigid-body motions in the axial direction. Theoretically, the Hw matrix given in Eqs. (14) will become non-invertible for a completely free shell. However, this numerical irregularity can be easily avoided by letting one of the bending-related springs have a very small stiffness, such as,k^4=106. Table 4 shows a comparison of the frequency parameters calculated using different techniques. While the results obtained from the current technique agree reasonably well with the other two reference sets, perhaps within the variance of different shell theories, the frequency parameters for the two lower order modes with rigid-body rotation (n=2 and 3) are clearly inaccurate which probably indicates that the inability of exactly satisfying the shell boundary conditions by the “beam functions” tends to have more serious consequence in such a case. Amazingly enough, the inextensional theory works very well in predicting the frequency parameters for the “rigid-body” modes (those with rigid-body motions in the axial direction). It is also seen that the frequency parameters of the rigid-body modes increases monotonically with the circumferential modal index n.

Mode n=2 n=3 n=4 n=5 n=6 n=7
(2.130) [22]
(2.132) [22]













Table 4.

Frequency parameters, Ω=ωRρ(1σ2)/E, for a free-free shell; R=0.5 m, L=4R, h=0.002R, and μ=0.28.

After it has been adequately illustrated how the classical boundary conditions can be easily and universally dealt with by simply changing the stiffness values of the restraining springs, we will direct our attention to shells with elastic end restraints. For the purpose of comparison, the problems previously studied in ref. [20] will be considered here. It was observed in that study that the tangential stiffness had the greatest effect on the natural frequency of the cylinder supported at both ends while the axial boundary stiffness had the greatest influence on the natural frequency of the cylinder supported at one end. It was also determined that natural frequencies varied rapidly with the boundary flexibility when the non-dimensionalized stiffness is between 10-2 and 102.

The frequency parameters for the “clamped”-free shell are shown in Table 5 for the reduced axial stiffness k^1L(1μ2)=1 (corresponding to ku*=1 in ref. [20]). It is seen that the current results are slightly larger than those taken from ref. [20]. The possible reasons include: 1) the difference in shell theories (the Flügge theory, rather than the Donnell-Mushtari, was used there), and 2) different Poisson ratios may have been used in the calculations.

Mode n=0 n=1 n=2 n=3 n=4 n=5



0.32866 (0.315*)

0.532604 (0.498)


Table 5.

Frequency parameters, Ω=ωRρ(1σ2)/E, for a “clamped”-free shell; R=0.00625 m, L=R, h=0.1R, μ=0.28, andk^1L(1μ2)=1.

Note: the numbers in parentheses are taken from ref. [20]

Although all eight sets of springs can be independently specified here, for simplicity we will only consider a simple configuration: a cantilevered shell with an elastic support being attached to its free (right) end in the radial direction. Listed in Table 6 are the four lowest natural frequencies for several different stiffness values. Obviously, the cases for k^7=0 and  represent the clamped-free and clamped-simply supported boundary conditions, respectively.

Mode k^7=0 k^7=0.01 k^7=0.1 k^7=1 k^7=1010
1319.99, ,

Table 6.

Natural frequencies in Hz for a clamped-elastically supported shell; L=0.502 m, R=0.0635 m, h=0.00163 m, E=2.1E+11, μ=0.28, ρ=7800 kg/m3;k^5=k^6=k^8=0.

The mode shapes for the three intermediate stiffness values are plotted in Figs. 2-4. It is seen that the modal parameters can be significantly modified by the stiffness of the restraining springs. The four modes in Fig. 2 for k^7=0.01 m-1 closely resemble their counterparts in the clamped-free case, even though the natural frequencies have been modified noticeably. While all the first four natural frequencies happen to increase, more or less, with the spring stiffness, the modal sequences are not necessarily the same. For example, when the spring stiffness k^7 is increased from 0.01 to 0.1 m-1, the third natural frequency goes from 886.66 to 926.17 Hz. However, this frequency drift may not necessarily reflect the direct effect of the stiffness change on the (original) third mode. It is evident from Figs. 2 and 3 that the third and fourth modes are actually switched in these two cases: the original third mode now becomes the fourth at 1200.88 Hz. It is also interesting to note that while stiffening the elastic support k^7 (from 0.01 to 0.1 m-1) has significantly raised the natural frequencies for the first two modes, the fourth mode is adversely affected: its frequency has actually dropped from 1023.61 to 926.17 Hz (see Figs. 2 and 3). A similar trend is also observed between the fourth mode for k^7=0.1 m-1 and the second mode for k^7=1 m-1, as shown in Figs. 3 and 4.

Figure 2.

First four modes for the clamped-elastically supported shell; k^7=0.01 m-1andk^5=k^6=k^8=0.

Figure 3.

First four modes for the clamped-elastically supported shell; k^7=0.1 m-1andk^5=k^6=k^8=0.

Figure 4.

First four modes for the clamped-elastically supported shell; k^7=1 m-1andk^5=k^6=k^8=0.

3.2. An exact solution based on the Flügge’s equations

To validate the exact solution method, the simply supported shell is considered again. Given in Table 7 are the calculated natural frequency parametersΩ=ωRρ(1σ2)/E. The current results agree well with the exact solutions based on Flügge’s theory [30], solutions based on beam functions [31] and three-dimensional linear elasticity solutions [30].

h/R n Ω=ωRρ(1σ2)/E
Ref. [30]a Ref. [31] Ref. [30]b Present
0.05 0 0.0929586 0.0929682 0.0929296 0.0929590
1 0.0161065 0.0161029 0.0161063 0.0161064
2 0.0393038 0.0392710 0.0392332 0.0393035
3 0.1098527 0.1098113 0.1094770 0.1098468
4 0.2103446 0.2102770 0.2090080 0.2103419
0.002 0 0.0929296 0.0929298 0.0929296 0.0929299
1 0.0161011 0.0161011 0.0161011 0.0161023
2 0.0054532 0.0054530 0.0054524 0.0054547
3 0.0050419 0.0050415 0.0050372 0.0050427
4 0.0085341 0.0085338 0.0085341 0.0085344

Table 7.

Comparison of values of the natural frequency parameter Ω=ωRρ(1σ2)/E for a circular cylindrical shell with simply supported boundary conditions, m = 1, R/l = 0.05, σ = 0.3.

a Exact solutions based on Flügge’s theory.

b Three-dimensional linear elasticity solutions.

The current solution method is also compared with the finite element model (ANSYS) for shells under free-free boundary condition. In the FEM model, the shell surface is divided into 8000 elements with 8080 nodes. The calculated natural frequencies are compared in Tables 8. An excellent agreement is observed between these two solution methods.

n m = 1 m = 2
FEM present difference (%) FEM present difference (%)
0 3229.8 3230.3 0.015% 5131.4 5131.1 0.006%
1 2478.6 2479.3 0.028% 4830.4 4830.6 0.004%
2 269.20 269.30 0.037% 276.62 278.58 0.704%
3 761.25 761.01 0.032% 770.99 771.62 0.082%
4 1459.2 1458.6 0.041% 1469.6 1469.3 0.020%
5 2359.4 2358.6 0.034% 2369.9 2369.0 0.038%

Table 8.

Comparison of values of the natural frequency for a circular cylindrical shell with free-free boundary conditions, L=0.502 m, R=0.0635 m, h=0.00163 m, σ =0.28, E=2.1E+11 N/m3, ρ=7800 kg/m3.

In most techniques, such as the wave approach, the beam functions for the analogous boundary conditions are often used to determine the axial modal wavenumbers. While such an approach is exact for a simply supported shell, and perhaps acceptable for slender thin shells, it may become problematic for shorter shells due to the increased coupling of the radial and two in-plane displacements. To illustrate this point, we consider relatively shorter and thicker shell (l=8R and R =39h). The calculated natural frequencies are compared in Table 9 for a clamped-clamped shell. It is seen that while the current and FEM results are in good agreement, the frequencies obtained from the wave approach (based on the use of beam functions) are significantly higher, especially for the lower order modes.

n m = 1 m = 2
FEM Ref. [32] present FEM Ref. [32] Present
0 3229.8 4845.5 3230.3 5146.0 8075.8 5139.8
1 1882.8 2350.2 1880.9 3850.7 4775.6 3848.9
2 899.59 985.48 898.18 2017.8 2303.4 2014.1
3 896.97 919.01 896.56 1390.9 1479.2 1388.9
4 1501.9 1517.45 1501.6 1676.4 1714.0 1676.0
5 2386.1 2402.05 2386.0 2472.5 2501.8 2472.6

Table 9.

Comparison of the natural frequencies for a circular cylindrical shell with clamped-clamped boundary conditions, L=0.502 m, R=0.0635 m, h=0.00163 m, σ=0.28, E=2.1E+11 N/m3, ρ=7800 kg/m3.

The exact solution method can be readily applied to shells with elastic boundary supports. Since the above examples are considered adequate in illustrating the reliability and accuracy of the current method, we will not elaborate further by presenting the results for elastically restrained shells. Instead, we will simply point out that the solution method based on Eqs. (27) is also valid for non-uniform or varying boundary restraint along the circumferential direction, which represents a significant advancement over many existing techniques.


4. Conclusion

An improved Fourier series solution method is described for vibration analysis of cylindrical shells with general elastic supports. This method can be easily and universally applied to a wide variety of boundary conditions including all the 136 classical homogeneous boundary conditions. The displacement functions are invariantly expressed as series expansions in terms of the complete set of trigonometric functions, which can mathematically ensure the accuracy and convergence of the present solution. From practical point of view, the change of boundary conditions here is as simple as varying a typical shell or material parameter (e.g., thickness or mass density), and does not involve any solution algorithm and procedure modifications to adapt to different boundary conditions. In addition, the proposed method does not require pre-determining any secondary data such as modal parameters for an “analogous” beam, or modifying the implementation algorithms to avoid the numerical instabilities resulting from computer round-off errors. It should be mentioned that the current method can be readily extended to shells with arbitrary non-uniform elastic restraints. The accuracy and reliability of the current solutions have been demonstrated through numerical examples involving various boundary conditions.



The authors gratefully acknowledge the financial support from the National Natural Science Foundation of China (No. 50979018).


  1. 1. 1964ForsbergK.Influence of boundary conditions on the modal characteristics of thin cylindrical shells,AIAA Journal ; 22150
  2. 2. ForsbergK.Axisymmetric and beam-type vibrations of thin cylindrical shells,AIAA Journal 19697221
  3. 3. WarburtonG. B.Vibrations of thin circular cylindrical shell, J. Mech. Eng. Sci. 19657399
  4. 4. WarburtonG. B.HiggsJ.Natural frequencies of thin cantilever cylindrical shellsJ. Sound Vib. 197011335
  5. 5. GoldmanR. L.Mode shapes and frequencies of clamped-clamped cylindrical shellsAIAA Journal1974121755
  6. 6. YuY. Y.Free vibrations of thin cylindrical shells having finite lengths with freely supported and clamped edges, J. Appl. Mech. 195522547
  7. 7. BerglundJ. W.KlosnerJ. M.Interaction of a ring-reinforced shell and a fluid mediumJ. Appl. Mech., 196835139
  8. 8. De SilvaC. N.TersteegC. E.Axisymmetric vibrations of thin elastic shellsJ. Acoust. Soc. Am. 19644666
  9. 9. SmithB. L.HaftE. E.Natural frequencies of clamped cylindrical shells,AIAA Journal19686720
  10. 10. LiX. A.newapproach.forfree.vibrationanalysis.ofthin.circularcylindrical.shellJ.Sound Vib. 20062969198
  11. 11. WangC.LaiJ. C. S.Prediction of natural frequencies of finite length circular cylindrical shellsApplied Acoustics200059385
  12. 12. ZhangX. M.LiuG. R.LamK. Y.Vibration analysis of thin cylindrical shells using wave propagation approach,Journal of Sound and Vibration 20012393397403
  13. 13. LiX. B.Study on free vibration analysis of circular cylindrical shells using wave propagationJournal of Sound and Vibration2008311667
  14. 14. ArnoldR. N.WarburtonG. B.Flexural vibrations of the walls of thin cylindrical shells having freely supported endsProc. Roy. Soc., A 1949197238
  15. 15. ArnoldR. N.WarburtonG. B.The flexural vibrations of thin cylindersProc. Inst. Mech. Engineers, A 19531676280
  16. 16. SharmaC. B.JohnsD. J.Vibration characteristics of a clamped-free and clamped-ring-stiffened circular cylindrical shell, J. Sound Vib. 1971;14459
  17. 17. SharmaC. B.JohnsD. J.Free vibration of cantilever circular cylindrical shells-a comparative study, J. Sound Vib. 197225433
  18. 18. SharmaC. B.Calculation of natural frequencies of fixed-free circular cylindrical shellsJ. Sound Vib. 19743555
  19. 19. SoedelW. A.newfrequency.formulafor.closedcircular.cylindricalshells.fora.largevariety.ofboundary.conditionsJ. Sound Vib. 198070309
  20. 20. LovedayP. W.RogersC. A.Free vibration of elastically supported thin cylinders including gyroscopic effects,J. Sound Vib, 1998217547562
  21. 21. AmabiliM.GarzieraR.Vibrations of circular cylindrical shells with nonuniform constraints, elastic bed and added mass: part I: empty and fluid-filled shellsJ. Fluids Struct, 200014669690
  22. 22. Leissa A. W. Vibration of Shells, Acoustical Society of America; 1993.
  23. 23. Qatu M.S.Recent research advances in the dynamic behavior of shells. part 2: homogeneous shells,Applied Mechanics Reviews200255415434
  24. 24. LiW. L.Free vibrations of beams with general boundary conditionsJ. Sound Vib. 2000237709
  25. 25. LiW. L.Vibration analysis of rectangular plates with general elastic boundary supportsJ. Sound Vib. 2004273619
  26. 26. LanczosC.Discourse on Fourier series.Hafner, New York; 1966
  27. 27. JonesW. B.HardyG.Accelerating convergence of trigonometric approximationsMath. Comp. 197024547
  28. 28. BaszenskiG.DelvosJ.TascheM. A.unitedapproach.toaccelerating.trigonometricexpansions.ComputMath. Appl. 19953033
  29. 29. BlevinsR. D.Formulas for Natural Frequency and Mode ShapeNew York: Van Nostrand Reinhold Company; 1979
  30. 30. KhdeirA. A.ReddyJ. N.Influence of edge conditions on the modal characteristics of cross-ply laminted shells, Computers and Structures 199034817
  31. 31. LamK. Y.LoyC. T.Effects of boundary conditions on frequencies of a multi-layered cylindrical shellJournal of Sound and Vibration1995363384
  32. 32. ZhangX. M.LiuG. R.LamK. Y.Vibration analysis of thin cylindrical shells using wave propagation approach,Journal of Sound and Vibration 20012393397403

Written By

Tiejun Yang, Wen L. Li and Lu Dai

Submitted: 07 March 2012 Published: 02 October 2012