Finite element for cylindrical thin shells under harmonic forces

Finite Elements in Analysis and Design 52 (2012) 83–92

Contents lists available at SciVerse ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

Finite element for cylindrical thin shells under harmonic forces Raydin Salahifar a,n, Magdi Mohareb b a b

´’ı´ Institute for Higher Education (BIHE), Tehran, Iran Department of Civil Engineering, Baha Department of Civil Engineering, University of Ottawa, Ottawa, ON, Canada K1N 6N5

a r t i c l e i n f o

abstract

Article history: Received 15 March 2010 Received in revised form 15 May 2011 Accepted 23 September 2011 Available online 7 December 2011

A super-convergent finite element is developed for the steady state analysis of circular cylinders under general harmonic forces based on thin shell theory. Simplifying assumptions in previous thin-shell formulations for straight cylinders are avoided. The resulting model is capable of capturing general in-plane and out-of-plain cross sectional distortions including ovalization, warping, radial extensibility, etc. A series of shape functions are then developed based on the exact solution of the field equations and used to formulate a cylindrical finite element capable of capturing thin shell behaviour with very few degrees of freedom. Through examples, the element is shown to efficiently model cylinders subject to in-phase and out-of phase harmonic loading and support excitation. & 2011 Elsevier B.V. All rights reserved.

Keywords: Finite element Circular cylindrical thin shell Pipe Harmonic forces Exact shape functions

1. Scope and literature review Circular cylindrical thin shells (CCTSs) are commonly used as structural members and as mechanical components. Also, they are key components of piping and conveyance systems. In these applications, CCTSs are commonly subjected to harmonic excitations induced by pumps, turbines, generators, etc. Under these loads, the steady-state response of the structure needs to be determined for fatigue design considerations since it represents the sustained effect. In contrast, the transient response of the structure is of little importance in fatigue design as it tends to dampen out quickly. In complex systems (e.g., long pipeline systems), shell based analytical solutions become impractical. Also, in such systems, finite element models based on conventional shell elements aimed at capturing localised shell effects necessitate a large number of degrees of freedom and are impractical as design tools. A promising class of finite element developments expresses the unknown displacement fields in terms of Fourier series with respect to the circumferential coordinate. In the longitudinal direction, various researchers have adopted different approximate shape functions such as Hermitian polynomials (e.g. [1–3]), Lagrangian interpolation functions (e.g. [3–6]), Harmonic functions (e.g. [2,7,8]), Chebyshev polynomials (e.g. [9]) or a combination thereof. In this class of finite elements, the mesh is discretized only along the longitudinal direction of the cylinder

and no discretization is needed in the circumferential direction, resulting in significant computational efficiency. The present formulation lies in this category of elements but departs from previous formulations in the following respects: (1) It adopts a family of interpolation functions for the displacement fields based on the closed form solution of the homogeneous form of the equilibrium equations. A similar approach was recently implemented for static loading [10] under the simplified shell theory in [11]. (2) While most of the solutions developed focus on either the static response (e.g. [10,12–14]) or on the free vibration (e.g., [15,16]), the present study extends the developments to the steady state response of CCTSs under harmonic loads. The solution directly captures the steady state response without the need to extract the modal shapes and natural frequencies. The formulation is capable of modelling multiple out-of-phase harmonic forces. (3) While most studies (e.g., [3,5,12,14]) typically use only a limited number of harmonic terms, the present formulation is applicable to any number of terms. In practicality, the number of terms in a series will be truncated when the accuracy sought is reached. The authors have provided in [17] a detailed tabular comparison of several studies in terms of number of harmonic terms used. (4) The proposed finite element is formulated based on a new tensor based thin shell theory derived by the authors in [18].

n

Corresponding author. Tel.: þ98 912 547 2292. E-mail addresses: [email protected], [email protected], [email protected] (R. Salahifar), [email protected] (M. Mohareb). 0168-874X/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2011.09.014

A comprehensive review of finite element formulations for pressure vessels and piping is provided in Mackerle [19–22].

84

R. Salahifar, M. Mohareb / Finite Elements in Analysis and Design 52 (2012) 83–92

The present formulation is aimed at developing a cylindrical element based on shell formulation. In this class of elements, the displacement fields are harmonic in the circumferential direction in order to enforce the continuity of displacements around the circumference. Hence, the following review focuses on shell based cylindrical elements and pipe elements, either curved or straight, with this particular feature. Various studies have used different interpolation schemes in the longitudinal direction and a limited number of harmonic modes as detailed in the following. Based on the principle of virtual work, Ohtsubo and Watanabe [1] developed a ring element for small strains for torus elements. Displacement components are interpolated by Hermitian polynomials along the elbow. Their model captures both ovalization and warping effects. Takeda et al. [2] modified the work of Ohtsubo and Watanabe [1] by including the transverse shear effects in their formulation and adopting 2D Lagrangian bilinear functions to interpolate the circumferential, and radial displacements and rotations. The longitudinal displacement is interpolated using Fourier series in the circumferential direction and linear functions in the longitudinal direction. Bathe and Almeida [4] developed a four node elbow finite element for in-plane flexure analysis. Based on Lagrangian cubic interpolation of the displacement fields, their model captures circumferential shear deformation and ovalization effects and neglects warping effects. They extended their model to investigate interaction effects between elbow elements and straight pipe segments [23]. Also, they added geometric non-linear effects to the formulation [24]. Militello and Huespe [6] extended the model by Bathe and Almeida [4] to capture warping deformations by adopting a Lagrangian cubic interpolation for the longitudinal displacement. The tangential and radial displacements were interpolated using harmonic functions. Abo-Elkhier [25] generalised the model by Bathe and Almeida [4] by adopting more complete displacement– strain relations into the formulation. In a ring formulation, Melo and Castro [26] used simple linear interpolation for displacements, which captures both in-plane and out-of-plane flexure of curved pipes. Sivadas and Ganesan [27] developed a finite element method for the vibration of a thin orthotropic cantilever circular cylindrical shell with variable thickness in axial direction. In their fivenode element, shear deformation effects are included, and fourth order Lagrangian polynomials are used to interpolate displacements along the element. Stanley and Ganesan [28] studied the response of cantilever CCTS with thickness discontinuities subject to harmonic axisymmetric loads including uniform internal pressure and circular ring load acting at cantilever tip. The longitudinal and circumferential displacements were linearly interpolated while the radial displacement was interpolated using cubic Hermitian functions. Xi et al. [29] developed a semi-analytical model for the vibration of composite shells of revolution based on the Reissner–Mindlin shell theory. The model is based on a two node straight conical frustum element capable of modelling shells of revolutions of general geometries. Yan et al. [30] developed an element for the prediction of the elasto-plastic limit load of pipe elbows by relaxing the pipe inextensibility assumption in the radial direction based on the model of Militello and Huespe [6]. They used cubic Lagrangian interpolation for the longitudinal displacement and Hermitian polynomials for the tangential and radial displacements. Pinto Correia et al. [3] developed a frustra conica finite element for the analysis of laminated axisymmetric shells based on a higher order shell theory. The radial displacements were interpolated using Hermitian polynomials while the longitudinal

and circumferential displacements were interpolated using linear Lagrangian interpolation. Karamanos [13] formulated an element that models buckling instability in pipes. Subsequent extension of the model involved adding plasticity and internal pressure effects (Karamanos et al. [14]). Santos et al. [5] developed a finite element model based on an axi-symmetric 3D formulation to study the bending, free vibration and buckling behaviour of shells of revolution made of laminated orthotropic elastic material. The displacement fields were linearly interpolated across the thickness and in the form of Lagrange polynomials in the longitudinal direction. Fonseca et al. [8] presented the development of two different finite piping elbow elements. The formulation is based on high order polynomial or trigonometric functions for rigid beam displacement and Fourier series for warping and ovalization phenomena of cross-tubular section.

2. Formulation 2.1. Geometry and loading The cylinder with mid-surface radius R, length L, thickness h and density r is subjected to body forces px, ps and pz respectively acting along the longitudinal, circumferential, and radial directions. Under the given forces, the cylinder is assumed to undergo mid-surface displacements u, v and w along the x, s and z coordinates, respectively. The sign convention adopted for the coordinates is illustrated in Fig. 1. 2.2. Assumptions Given the frequent to permanent nature of harmonic loads in the aforementioned applications, CCTSs are normally designed to respond to harmonic effects within the linearly elastic range of their material. Thus, the cylinder material is assumed to be linearly elastic. The formulation is restricted to cylindrical shells of circular cross-sections. The Love–Kirchhoff kinematic assumptions for thin shells are adopted, i.e, a straight line originally normal to the mid-surface remains normal. Displacements and strains are assumed small and stresses normal to the shell midsurface and transverse shear stress are neglected. The formulation is general in that it captures all membrane and bending stresses, cross-sectional distortions (such as ovalization) and warping effects under general loads and boundary conditions.

Fig. 1. Cylinder geometry and coordinate definition.

R. Salahifar, M. Mohareb / Finite Elements in Analysis and Design 52 (2012) 83–92

2.3. Expressions for applied forces Harmonic body forces applied to a cylindrical shell are expressed as þ1 X

px ðx,s,z,tÞ ¼ px ðx,s,zÞeiot ¼

pxn ðx,zÞeins=R eiot

n ¼ 1 þ1 X

ps ðx,s,z,tÞ ¼ ps ðx,s,zÞeiot ¼

þ1 X

2.5. Variational expression and equilibrium conditions pzn ðx,zÞeins=R eiot

ð1a2cÞ

n ¼ 1

in which, o is the exciting angular frequency. In Eqs. (1a–c), the continuity of the applied body forces in the circumferential direction is ensured by assuming the dependency of px, ps, pz on the circumferential coordinate s to take the harmonic form eins/R. Other load distributions such as surface tractions, ring loads, line loads or point loads can always be conceived as special cases of Eqs. (1a-c) using the Dirac Delta function as needed. When the acting forces are in phase, force distributions px ðx,s,zÞ, ps ðx,s,zÞ, pz ðx,s,zÞ are functions with real values. However, the current model is capable of taking complex force distributions, in order to model forces with the same angular frequency o but with different phase angles [18].

vðx,s,tÞ ¼ vðx,sÞe

¼

þ1 X

! Bn ðxÞe

n ¼ 1

wðx,s,tÞ ¼ wðx,sÞeiot ¼

ins=R

þ1 X

d

Y

Z

dðT n V n Þdt þ

t1

Z

t2

dW n dt ¼ 0

where P is the mechanical energy of the system, Tn is the kinetic energy, Vn is the internal strain energy and Wn is the external work on the system. Symbol d denotes the variation operator and integration is performed over an arbitrary period starting from time t1 to time t2. For thin shells, variations dTn, dVn and dWn are ([31]) Z

dW n ¼

i

rU_ dU_ i dV, dV n ¼

Z V

Z V

Eabld Zld dZab dV,

F i dU i dV

ð4a2cÞ

in which indices i,j¼1,2,3, indices a,b,l,d ¼1,2, the differential pffiffi volume element dV ¼ Jdxdsdz and J denotes the determinant of the Jacobian. Dots on top of the symbols denote differentiation with respect to time and integrations are performed over the shell volume. Symbols Ui and U i ¼ g ij U j respectively denote the covariant and contravariant components of the displacement fields, and gij are the components of the contravariant metric tensor. The physical 1

3

U 3 ¼ U ¼ W. Symbol Eabld is the contravariant constitutive elasticity tensor and Zld is the covariant strain tensor. The non-zero 1111

ð2a2cÞ

n ¼ 1

In Eqs. (2a–c), the use of the harmonic continuous functions eins/R ensures that the displacement fields and their derivatives with respect to the circumferential coordinate are continuous. The displacement expressions in Eqs. (2a–c) are quite general. The terms corresponding to n¼0 represent the axi-symmetric components of the displacements (i.e., A0(x) represents the longitudinal displacement characteristic of a cylinder under longitudinal end traction uniformly distributed around the circumference, B0(x) is the tangential displacement characteristic of a cylinder under St. Venant torsion, and C0(x) is the radial extensibility characteristic of a cylinder under internal or external pressure). The terms n¼ 71 represent beam-like displacements: (a) Combinations of displacements v and w based on n ¼ 71 were observed to yield a nearly circular deformed cross-sections, and (b) the deformed configuration of the cylinder cross-section undergoing the longitudinal displacement u based on n ¼ 71 remain co-planar, in line with conventional beam theories.

2

components of Ui and Ui are U 1 ¼ U ¼ U, U 2 ¼ U ¼ V and

2211

!

ð3Þ

t1

physical components of Eabld are E

eiot

C n ðxÞeins=R eiot

t2

¼

V

In line with the assumptions in Section 2.2 the longitudinal, circumferential and radial displacements U, V and W of an arbitrary point within the shell are expressed in terms of the mid-surface displacements u, v and w through U¼u zw,1, V¼v(1þ(z/R))  zw,2 and W¼w, in which subscripts ‘‘,1’’ and ‘‘,2’’ denote differentiation with respect to x and s, respectively. The present study is aimed at capturing only the steady state response of the system. Under the given harmonic forces, the mid-surface displacements corresponding to the steady state response are also harmonic, i.e., ! þ1 X iot ins=R uðx,s,tÞ ¼ uðx,sÞe ¼ An ðxÞe eiot n ¼ 1

In general, the variational form of Hamilton’s principle is expressed as

dT n ¼

2.4. Steady state displacement functions

iot

The terms 9n9Z2 represent cross-sectional distortions (specifically, displacements v and w based on n ¼ 72 represent the cross-section ovalization, while higher order terms represent higher modes of cross-section distortions). Also, displacement u with 9n9 Z2 represent the deviations from the co-planar behaviour represented by mode n ¼ 71 and thus represent warping deformations.

psn ðx,zÞeins=R eiot

n ¼ 1

pz ðx,s,z,tÞ ¼ pz ðx,s,zÞeiot ¼

85

2

1212

2222

¼E

1221

1122

¼ E=ð1n2 Þ, E

2112

¼

2121

E ¼ Ev=ð1n Þ and E ¼E ¼E ¼E ¼ E=2ð1þ nÞ where E is Young’s modulus and n is Poisson’s ratio. The physical components of Zld are Z11 ¼ u,1 zw,11 , Z22 ¼ v,2 þðw=ðR þzÞÞ ðzw,22 =ð1 þðz=RÞÞ and Z12 ¼ Z21 ¼ ½ðu,2 Þ=ð1 þðz=RÞÞþ v,1 ð1 þðz=RÞÞ zw,12 ðð2R þzÞ=ðR þ zÞÞ=2. In Eq. (4c), Fi are the contravariant compo1

2

nents of the external body forces, where F ¼ px , F ¼ ð1 þðz þ RÞÞps 3

and F ¼ pz are their physical components. After performing integration by parts with respect to t, x and s, integrating across the thickness of the shell (along the z coordinate) and around the circumference (along the s coordinate), and evoking the stationary condition of the Hamilton’s functional, it can be shown ([18],[31]) that

d

Y

Y Y Y Y Y ¼ d IF þ d BC þ d BF þ d FE þ d MF ¼ 0

ð5Þ

in which, the time boundary term dPIF ¼0 provides the initial and final conditions of the system at arbitrary time limits t1 and t2. This term is not relevant to steady-state response of the system. The terms dPBC þ dPBF provide the boundary condition equations of the system at x ¼0 and x ¼L, where the term dPBF is contributed by member forces. The terms dPFE þ dPMF represent the field equations of the system where the term dPMF is the

86

R. Salahifar, M. Mohareb / Finite Elements in Analysis and Design 52 (2012) 83–92

contribution of the member forces. In Eq. (5), the following terms are defined

satisfy the homogeneous part of the field equations (Eq. 6), i.e.

ð8Þ It can be shown ([18,31]) that the solution of Eq. (8) is An ðxÞ ¼

8 X

A~ nj Anj emnj x ,

Bn ðxÞ ¼

j¼1

C n ðxÞ ¼

8 X

A~ nj Bnj emnj x ,

j¼1

8 X

A~ nj C nj emnj x

ð9a2cÞ

j¼1

ð6a2dÞ in which Z þ ðh=2Þ   z z xw pxn 1 þ dz, f n ¼  pxn z 1 þ dz R R ðh=2Þ ðh=2Þ Z þ ðh=2Þ  Z þ ðh=2Þ  z 2 z sv sw s s fn ¼ pn 1þ dz, f n ¼  pn z 1 þ dz R R ðh=2Þ ðh=2Þ Z þ ðh=2Þ  z zw z ð7a2eÞ pn 1 þ dz fn ¼ R ðh=2Þ xu

fn ¼

Z

þ ðh=2Þ

In Eq. (6), symbol D is the differential operator and Dd (d¼2,3,4) denotes the dth differentiation of the argument function. Constants ij H, ij Hn , ij G and ij Gn are defined in Appendix A. As a convention, the left superscript i (i ¼1,y,4) of constants ij H, i i i j H n , j G and j Gn indicate the constants are related to the ith boundary or integral term, as applicable. The left subscript j denotes the jth coefficient of ith boundary or integral term. The presence of subscript n in constants ij Hn and ij Gn denotes that the constants are dependent on the circumferential mode number n. This alternative form of the variational form of the Hamilton’s principle (Eq. 5) has the following advantages: (a) It provides a unified and consistent system of field equations and corresponding boundary condition equations. (b) With the methodology employed in the current study, as it will be discussed in Section. 2.6, it yields a more concise system of equations for the discretized formulation. It is noted that the present field equations and corresponding boundary conditions based on tensor calculus (Eq. (6)) differ from previous theories (e.g., Flugge and Donnell–Mushtari–Vlasov (DMV) [32], Timoshenko [11], Niordson and Morley-Koiter [33], Novozhilov [34], Saada [35], etc.), in that they avoid simplifying assumptions (such as 1þ z=R  1) beyond those cited in Section 2.2 of the paper. As a result, the present shell formulation yields more general field equations and boundary conditions. It can be shown [31] that the field equations of Flugge and DMV can be recovered as special cases of the present field equations. 2.6. Formulating the shape functions In Eq. (2), functions An(x), Bn(x) and Cn(x) are unknown. Conventionally, approximate polynomial interpolation functions would be used at this stage to relate functions An(x), Bn(x) and Cn(x) to nodal displacements. The present study departs from this convention by adopting interpolation functions, which exactly

in which subscript n is the mode number in Fourier expansion series in circumferential direction (Eq. (2)). Integration constants A~ nj (j¼1,y,8) are unknowns to be determined from boundary conditions and mnj are the eight roots of the characteristic equation an m8n þ bn m6n þcn m4n þ dn m2n þ en ¼ 0 in which coefficients an,y,en are defined in Appendix B. It is observed that the characteristic roots of modes n and  n are identical. Constants Anj , Bnj and C nj are defined as 8 > < 1 ; n ¼ 0, j ¼ 1,. . .,6 Anj ¼ 0 ; n ¼ 0, j ¼ 7,8 ð10Þ > : 1 ; n a0, j ¼ 1,. . .,8 8 0 > > < Bnj ¼ 1 > > K nj:13 K nj:21 K nj:11 K nj:23 : K nj:12 K nj:23 K nj:13 K nj:22

8 1 > > < C nj ¼ 0 > K K K K > : nj:11 nj:22 nj:12 nj:21 K nj:12 K nj:23 K nj:13 K nj:22

; n ¼ 0, j ¼ 1,. . .,6 ; n ¼ 0, j ¼ 7,8

ð11Þ

; n a0, j ¼ 1,. . .,8 ; n ¼ 0, j ¼ 1,. . .,6 ; n ¼ 0, j ¼ 7,8

ð12Þ

; n a 0, j ¼ 1,. . .,8

in which 1 1 K nj 1 3 K nj

¼ 11 Gn þ 12 Gm2nj ,

2 3 K nj

¼ 32 K nj ¼ 24 Gn þ 25 Gn m2nj ,

¼

3 1 K nj

1 1 2 2 K nj ¼ 1 K n,j ¼ 3 Gn mnj 1 1 2 3 ¼ 4 Gn mnj þ 5 Gmnj , 2 K nj ¼ 22 Gn þ 23 Gm2nj 3 3 K nj

¼ 35 Gn þ 36 Gn m2nj þ 37 Gm4nj

ð13a2fÞ

In a matrix form, Eqs. (9a–c) are expressed as An ðxÞ ¼ /SAn ðxÞS18 fAn g81 Bn ðxÞ ¼ /SBn ðxÞS18 fAn g81 C n ðxÞ ¼ /SCn ðxÞS18 fAn g81

ð14a2cÞ

in which

ð15a2cÞ and ð16Þ In a finite element context, all unknown coefficients {An}8  1 are to be expressed in terms of nodal displacements An(0), Bn(0), Cn(0), C 0n ð0Þ, An(L), Bn(L), Cn(L) and C 0n ðLÞ. From Eqs. (14a–c), one has fDn g81 ¼ ½Ln 88 fAn g81

ð17Þ

in which

ð18Þ

R. Salahifar, M. Mohareb / Finite Elements in Analysis and Design 52 (2012) 83–92

and ð19Þ where

ð20Þ Solving Eq. (17) for the unknown coefficients {An}8  1 yields fAn g81 ¼ ½Ln 1 88 fDn g81

ð21Þ

87

analytical solution is developed for static loads, were solved using the present formulation and the results were found in excellent agreement with those in [10]. The maximum difference was less than 2%. As will be explained under Example 3, convergence to static solutions was attained by applying an exciting frequency significantly smaller than the fundamental frequency of the structure. Also, all examples in [18], in which a shell based analytical solution for harmonic forces is developed, have been solved with the proposed finite element and all results are found to exactly coincide with those reported in [18].

By substituting Eq. (21) into Eq. (15), one obtains An ðxÞ ¼ /NAn ðxÞS18 fDn g81 Bn ðxÞ ¼ /NBn ðxÞS18 fDn g81 C n ðxÞ ¼ /NCn ðxÞS18 fDn g81

4. Examples ð22a  cÞ

in which /NAn ðxÞS18 ¼ /SAn ðxÞS18 ½Ln 1 88 /NBn ðxÞS18 ¼ /SBn ðxÞS18 ½Ln 1 88 /NCn ðxÞS18 ¼ /SCn ðxÞS18 ½Ln 1 88

ð23a  cÞ

It is noted that /NAn(x)S1  8, /NBn(x)S1  8 and /NCn(x)S1  8 are the exact shape function row vectors. In a cylinder consisting of multiple finite elements, the use of the shape functions in Eqs. (23a–c) while satisfying the interelement continuity conditions (arising from meeting the dPBC þ dPBF ¼0 in Eq. (5)) ensures C0 continuity in the longitudinal direction for functions An(x),Bn(x), and C1 continuity for function Cn(x). From Eqs. (2a–c), one recalls the identities ðu,v,wÞ ¼

þ1 X

½An ðxÞ,Bn ðxÞ,C n ðxÞeins=R

n ¼ 1

At an inter-element interface, given that the same Fourier expansion term eins/R is used to characterize the displacement contributions from mode n, C0 inter-element continuity for the longitudinal and tangential displacements, u and v, and C1 continuity for the radial displacement w are ensured for any generator along the cylinder.

4.1. Example 1 A 5 m long cantilever beam composed of two cylindrical segments of radius R¼0.2 m is subject to a harmonic point load Pq sin ot assumed to act on the middle surface (Fig. 2), with an exciting frequency o ¼ 10  2p rad=s and Pq ¼1000 N. Dimensions are given on Fig. 2. It is required to calculate the midsurface displacements of the generator of the beam at a ¼301 (i.e., s/R¼11p/6). For the present and following examples, material is steel with E¼200 GPa and n ¼0.30. Considering the coordinate configuration defined in Fig. 1, the time-independent contributions of the applied body forces are px ðx,s,zÞ¼0, ps ðx,s,zÞ¼P q sin a Diracðx5ÞDiracðsð11pR=6ÞÞ DiracðzÞ and pz ðx,s,zÞ ¼ P q cos a Diracðx5ÞDiracðsð11pR=6ÞÞ DiracðzÞ. The nodal displacements An(0), Bn(0), Cn(0) and C 0n ð0Þ at the clamped end (x¼0) of the beam are set to zero for all modes in order to model the fixity conditions. One element is used for each segment in order to model the change in thickness. The problem is solved twice: for 7 rnr þ7 and  10rn r þ10. The mid-surface displacements of the generator at s/R¼11p/6 are plotted in Fig. 3 and compared to the results based on the Abaqus

2.7. Discretized equilibrium equations From Eqs. (22a–c), by substituting into Eqs. (6a–d) and then Eq. (5), and noting that dPIF ¼0 and the shape functions derived exactly satisfy the homogeneous part of the field equations (i.e. the term dPFE vanishes when the exact shape functions are used, in contrast to other shape function selections, which do not vanish dPFE), one obtains þ1 X

dfDn gT18 ð½Kn 88 fDn g81 fFn g81 Þ ¼ 0

Fig. 2. Example 1, Cantilever cylinder with tip point load.

ð24Þ

n ¼ 1

in which the entries of the stiffness matrix contribution [Kn]8  8 arise from substituting Eqs. (22a–c) into the expression of dPBC (Eq. 6a) and the energy equivalent load vector contribution {Fn}8  1 arise from substituting Eqs. (22a–c) into the expression dPBF þ dPMF (Eqs. (6b,d)). Both [Kn]8  8 and {Fn}8  1 correspond to the nth mode contributions and are defined in Appendix C. It is noted that the use of the shape functions, which exactly meet the homogeneous equilibrium conditions, eliminate discretization errors typically arising in finite element formulations.

3. Verifications against previous analytical solutions In order to assess the validity of the proposed finite element formulation, all examples reported in [10], in which a shell based

Fig. 3. Longitudinal, circumferential and radial displacements of the generator of the beam at s/R ¼11p/6.

88

R. Salahifar, M. Mohareb / Finite Elements in Analysis and Design 52 (2012) 83–92

model. In the Abaqus model two different meshes are used, a coarser mesh consisting of 72  287 ¼20,664 S4R shell elements and a finer mesh consisting of 144  573¼82,512 elements. Fig. 3 shows that all four solutions are indistinguishable. The model developed in the present study is highly accurate and efficient with far fewer degrees of freedom (two elements with a total of 64 DOFs for  7rn r þ7 in the current study compared to 20,664 elements with 123,984 DOFs in the Abaqus coarser mesh). A 3D plot of the deformed configuration of the cylinder based on the current model is shown in Fig. 4a, which shows the ability of the proposed model to capture localised effects. It is noted that the fine mesh depicted in Fig. 4a is provided for a clear visualisation of the deformed configuration and is no indicator of the number of elements adopted in the solution based on the present formulation. A plot based on the finer Abaqus mesh is also shown on Fig. 4b for comparison.

4.2. Example 2 A steel pipe (Fig. 5) of radius R¼0.25 m and thickness h¼5 mm is fully restrained at its left support x¼0. At the right end (x¼5 m) the pipe cross-section is connected to a rigid shear diaphragm, which prevents in-plane cross-sectional distortions and out-of-plane warping. The right end of the pipe is free to undergo transverse displacements but is restrained in the lateral and longitudinal directions. It is required to determine the response of the system,

Fig. 5. Pipe configuration.

if the right support is subject to a vertical harmonic excitation of frequency o ¼ 20  2p rad=s and a 10 mm amplitude. At the clamped end of the pipe (x¼ 0), all nodal displacements An(0), Bn(0), Cn(0) and C 0n ð0Þ vanish. At the excited support (x¼5), the given harmonic support excitation can be mathematically described as vðx ¼ 5,s,tÞ ¼ ½0:01 sinðs=RÞeiot and wðx ¼ 5,s,tÞ ¼ ½0:01 cosðs=RÞeiot . Considering Eqs. (2a–c), these exciting displacements contribute only to modes n ¼ 71. Thus, the nodal displacements are A1 ð5Þ ¼ A þ 1 ð5Þ ¼ C 01 ð5Þ ¼ C 0þ 1 ð5Þ ¼ 0 and B  1(5)¼B þ 1(5) ¼  C  1(5)¼  C þ 1(5)¼0.005 mm. Since all nodal displacements An(0), Bn(0), Cn(0), C 0n ð0Þ, An(5), Bn(5), Cn(5), and C 0n ð5Þ vanish for modes na 71, Eqs. (22a–c), which provide the displacement fields for the case of zero member forces imply that the displacement fields An(x),Bn(x),Cn(x) vanish and the only non-zero displacement contributions arise from modes n ¼ 71. The longitudinal and circumferential stresses at the top three generators (inner surface z¼  h/2, mid-surface z¼0 and outer surface z¼ þh/2) of the beam are plotted in Fig. 6a and b, and the maximum shear stresses at the beam mid-height plotted in Fig. 6c. In all cases, comparisons are made to Abaqus results. The Abaqus model consists of 71,520 S4R shell elements with a finer mesh around the supports. The difference between both sets of results is indistinguishable within the scale of the figures. It is noted that the current model has reached this accuracy by using only eight degrees of freedom. 4.3. Example 3 It is required to obtain the static response of a fixed-fixed horizontal pipe, with the same dimensions and material properties as in Example 2, to its self weight. The present formulation yields a singular stiffness matrix for a vanishing exciting frequency. However, the problem can be approximately solved by considering a very low exciting frequency of o ¼ 2p rad=s (which is less than 1% of the first natural frequency of the pipe) in order to simulate static loading. The mid-span vertical deflection at the mid-surface for a point lying on top generator was found to be 0.0260 mm, which is identical to that based on a static analysis in Abaqus with 81,600 S4R shell element. The excellent agreement obtained is an indication of the ability of the formulation to capture the static response as a special case. 4.4. Example 4

Fig. 4. Deformed Configuration for a cylinder subject to point load based on (a) current Model (solid lines denote the top and bottom generators and slashdotted line denotes the generator at s/R ¼11p/6), and (b) Abaqus S4R shell model.

A horizontal cylindrical thin shell of span L¼ 6 m, radius R¼0.1 m and thickness h¼0.003 m is fixed at both ends. The shell is subject to a distributed harmonic patch load acting at mid-span within a rectangular area bounded with the lines s/R¼ 7 p/8 and x¼370.06 (Fig. 7a). Two types of uniform harmonic loads act on the patch: (a) A vertical traction load qV ¼ 1,000 eiot KN=m2 and (b) a longitudinal load qL ¼ 1,000 eiðot þ y1 Þ KN=m2 , which is y1 out-of-phase relative to qV. Such out-of-phase loading may be induced by rotary equipment with an eccentric

R. Salahifar, M. Mohareb / Finite Elements in Analysis and Design 52 (2012) 83–92

89

Fig. 6. Maximum stresses along the top generator for inner surface (z ¼  h/2), mid-surface (z ¼ 0) and outer surface (z ¼ þh/2). Solid lines denote results based on the current study and dotted lines denote results based on the Abaqus.

The mathematical expressions of the corresponding body forces are px ðx,s,z,tÞ ¼ P q ðx,s,zÞqL eiðot þ y1 Þ ps ðx,s,z,tÞ ¼ Pq ðx,s,zÞ½qV sinðs=RÞeiot pz ðx,s,z,tÞ ¼ P q ðx,s,zÞ½qV cosðs=RÞeiot

Fig. 7. (a) Patch loading area and the direction of the out of phase applied tractions qL and qV, and (b) distribution of the vertical traction qV.

mass. It is required to determine the steady state displacements under the combined effect of the two given loads with y1 ¼ p/4 rad. The exciting frequency is o ¼ 2p  60 rad=s.

ð25Þ

in which, Pq(x,s,z) ¼[H(x  2.94)  H(x  3.06)][1  H(s  pR/8)þH (s  15pR/8)]d(z  h/2) and H is the Heaviside unit step function. It is noted that circumferential and radial components of the applied forces result from decomposing the vertical traction force (Fig. 7b) into the radial and tangential directions. The timeindependent parts of the body forces are px ðx,s,zÞ ¼ qL Pq ðx,s,zÞ eiy1 , ps ðx,s,zÞ ¼ ½qV sinðs=RÞPq ðx,s,zÞ and pz ðx,s,zÞ ¼ ½qV cosðs=RÞ Pq ðx,s,zÞ. The number of Fourier terms needed to expand the loads in the circumferential direction was numerically found to be 15 (  7rn r þ7). The real components of the mid-surface displacement components along the top generator are plotted and compared to the results based on Abaqus in Fig. 8. The Abaqus shell model consisted of 96  916 ¼87,936 S4R shell elements. In the solution based on the current study, two meshes are used: (a) A single element along the patch load area and an element on

90

R. Salahifar, M. Mohareb / Finite Elements in Analysis and Design 52 (2012) 83–92

Appendix A. Coefficients of boundary and field terms The coefficients of the boundary and pffiffiffiffiffiffiffifield terms in Eqs. (6) are defined below. We recall that i ¼ 1, Fðh=RÞ ¼ ln½ð2 þ ðh=RÞÞ= ð2ðh=RÞÞ=ðh=RÞ and n is the circumferential mode number. 1 1H

¼ 1,

2 1 Hn

¼

1 2 Hn

¼  Rn in,

1n 1 in, 2 R

1 3H

1 4H

¼  Rn ,

2

2 2H

¼

h2 12R

¼ !

2

1n h 1þ 2 , 2 4R

2 3 Hn

1n h in 2 4R2

¼

2

2

1 1n 2 1n2 h h n ro2 , 32 H ¼  ¼  ðF1Þ R 2 E 12R 12R  2  h 3n 3 in H ¼  3 n 2 12R2 ( ) 2 2 h 3n 1n 2 1n2 h 3 þ H ¼  ð F 1 Þ n þ ro2 4 n 2 2 2 E 12 12R 3 1 Hn

2

h 1n2 , 36 H ¼ 12 Eh 2 2 h h 4 4 , 2 Hn ¼ nin, 1H ¼ 12R 12R2 3 5H

Fig. 8. Maximum mid-surface displacement components along the top generator.

¼

2

4 3 Hn

2

h

¼

12R

4 4H

nn2 ,

2

¼

Table 1 Mid-surface displacement components of Point J (in mm).

Current study Abaqus % Difference

Longitudinal

Circumferential

Radial

1 1 Gn

0.0600e0.857i 0.0586e0.859i 2.3

0.629e3.127i 0.621e3.128i 1.3

1.781e6.275i 1.754e6.275i 1.5

1 4 Gn

each side of the patch load. (b) Two elements along the patch load area and an element on each side of the patch load. The components of displacement of point J on the mid-surface at the corner of the patch load (Fig. 7a) are given and compared to those based on the Abaqus in Table 1. The magnitude of the imaginary exponent represents the out-of-phase angle with the applied tractions ps and pz. It is of interest to note that in the present example and all previous ones, the computational time based on the present study is three orders of magnitude faster than the solutions provided by Abaqus, while attaining the same degree of accuracy. 5. Conclusion A family of shape functions is formulated based on the exact solution of the differential equations of equilibrium of circular cylindrical thin shells subjected to harmonic forces. A super-convergent cylindrical finite element is developed based on the formulated shape functions. The element is capable of capturing warping, ovalization, radial extensibility, etc. The element efficiently captures the response of cylindrical elements and pipe systems under general harmonic forces, boundary conditions, and support excitations with very few degrees of freedom. Through four examples, the model is shown to predict the displacements and stresses in excellent agreement with those based on established shell finite element solutions. The computational efficiency of the element makes it particularly suitable for the analysis of complex systems such as long pipeline systems subject to harmonic forces.

¼

1

F

2

R

1n 2 1n2 n þ ro2 , 2 E

1 2G

1 3 Gn

¼ 1,

1 1þn in R 2

2

1 1n 2 1n2 h  ðF1Þ n  , ro2 R R 2 E 12R 2 1 n 1 6G ¼ Eh ¼

¼

n

2

1 5G

¼

h 12R

2

1 1þn 1 1n2 h in, 22 Gn ¼  2 n2 þ ro2 1 þ R 2 E 12R2 R ! 2 2 h 1n 1 1n2 h 2 , 24 Gn ¼ 2 in G ¼ 1 þ ro2 2 in 3 2 2 E R 4R 6R 2 1 Gn

¼

2 5 Gn

¼

2

3n in, 12R2 2 h

2 6G

¼

h 12 ðA:1a2qÞ

!

1n2 Eh 2

2

1 1n 2 1n2 h h ðF1Þ n þ , 32 G ¼ ro2 R R 2 E 12R 12R 2 2 1 1n2 h h 3n 3 in ro2 2 in, 34 Gn ¼ 3 Gn ¼  2 in þ E R 6R 12R2 2 1 2 1 3 2 4 5 Gn ¼  2 F þ 2 ðF1Þn  2 ðF1Þn R R R 2 1n2 1n2 h þ ro2 þ ro2 n2 2 E E 12R " # 2 2 h 3þn 1n 2 1n2 h 3 þ ð G ¼ F 1Þ ro2 n  6 n 2 2 2 E 12 12R 3 1 Gn

¼

n

þ

2

3 7G

¼

h , 12

3 8G

¼

1n2 , Eh

3 9 Gn

¼

1n2 in, EhR

3 10 G

¼

1n2 Eh ðA:2a2vÞ

Also, 2 1 Gn

¼ 13 Gn ,

3 1 Gn

¼ 14 Gn ,

3 2G

¼ 15 G,

3 3 Gn

¼ 24 Gn ,

3 4 Gn

¼ 25 Gn ðA:3Þ

Appendix B. Constants in characteristic equations Acknowledgement The authors would like to express their deepest gratitude to the Baha´’ı´ Institute for Higher Education (BIHE) for its encouragement and support to the first author. Also, the research grant support from the Natural Science and Engineering Research Council of Canada (NSERC) to the second author is gratefully acknowledged.

The coefficients in characteristic equation of the system discussed in Section 2.6 are defined below. an ¼ ð11 Gn Þð22 Gn Þð35 Gn Þð11 Gn Þð24 Gn Þ2 bn ¼ ½ð35 Gn Þð13 Gn Þ2 þ2ð13 Gn Þð14 Gn Þð24 Gn Þð22 Gn Þð14 Gn Þ2 ð14 GÞð24 Gn Þ2 2ð11 Gn Þð25 Gn Þð24 Gn Þ þð12 GÞð22 Gn Þð35 Gn Þ þð23 Gn Þð11 Gn Þð35 Gn Þ

R. Salahifar, M. Mohareb / Finite Elements in Analysis and Design 52 (2012) 83–92

þ ð11 Gn Þð22 Gn Þð36 Gn Þ

91

Force vector

cn ¼ ½ð36 Gn Þð13 Gn Þ2 þ 2ð13 Gn Þð14 Gn Þð25 Gn Þ þ2ð15 GÞð24 Gn Þð13 Gn Þ

The force vector results from Eqs. (6b,d) (i.e. dPBF and dPMF ). The term dPBF can be expressed as

ð23 Gn Þð14 Gn Þ2 2ð15 Gn Þð22 Gn Þð14 Gn Þð11 Gn Þð25 Gn Þ2 2ð12 Gn Þð24 Gn Þð25 Gn Þ þ ð12 GÞð23 GÞð35 Gn Þ

xw

xw

½36 H dC n ðxÞf n ðxÞL0 ¼ ½36 HðdfDn gT18 /NCn ðxÞST81 Þf n ðxÞL0

þ ð12 GÞð22 Gn Þð36 Gn Þ þ ð23 GÞð11 Gn Þð36 Gn Þ þ ð37 GÞð11 Gn Þð22 Gn Þ

xw ¼ ½dfDn gT18 ð36 H/NCn ðxÞST81 f n ðxÞÞL0 3 T ¼ dfDn g18 f6 f n g81

dn ¼ ½ð22 Gn Þð15 GÞ2 þ 2ð15 GÞð13 Gn Þð25 Gn Þ2ð23 GÞð14 Gn Þð15 GÞð37 GÞð13 Gn Þ2 ð12 Gn Þð25 Gn Þ2 þð12 Gn Þð23 GÞð36 Gn Þ þð12 GÞð37 GÞð22 Gn Þ þ ð23 GÞð37 GÞð11 Gn Þ en ¼ ð15 GÞ2 ð23 GÞ þð12 GÞð23 GÞð37 GÞ

ðB:1a2eÞ

Appendix C. Entries of stiffness matrix and force vector

where 3

Z

Stiffness matrix

xw

f6 f n g81 ¼ ½36 H/NCn ðxÞST81 f n ðxÞ10

0

L

ðE:3Þ

The first integral term in the first row of dPMF is Z L xu xu T 1 1 T 6 G dAn ðxÞf n ðxÞdx ¼ 6 GðdfDn g18 /NAn ðxÞS81 Þf n ðxÞdx 0

1

The entries of the stiffness matrix contribution corresponding to mode n result from the expansion of the nth terms of dPBC in Eq. (5). For each mode n, there are 16 contributions resulting from the matrix multiplication

¼ dfDn gT18 f6 f n g81 where 1

f6 f n g81 ¼ 16 G

Z

L 0

xu

/NAn ðxÞST81 f n ðxÞdx

In the same manner other similar terms yield Z L 2 sv /NBn ðxÞST81 f n ðxÞdx f6 f n g81 ¼ 26 G 0 Z L 3 xw0 f8 f n g81 ¼ 38 Gn /NCn ðxÞST81 f n ðxÞdx 0

3 f9 f n g81

The first of the 16 terms is ½ 11 H dAn ðxÞA0n ðxÞL0

¼

h

¼ 39 Gn

3

T 1 0 T 1 HðdfDn g18 /NAn ðxÞS81 Þð/NAn ðxÞS18 fDn g81 Þ

iL 0

¼ dfDn gT18 ½ 11 H/NAn ðxÞST81 /N0An ðxÞS18 L0 fDn g81

Z

f10 f n g81 ¼ 310 G

L

0

Z

0

sw

/NCn ðxÞST81 f n ðxÞdx

L

zw

/NCn ðxÞST81 f n ðxÞdx

ðE:4a2eÞ

The summation of Eqs. (E.3) and (E.4a–e) yields the force vector

1

¼ dfDn gT18 ½ 1 kn 88 fDn g81

fFn g81 ¼

where 1

½1 kn 88 ¼ ½11 H/NAn ðxÞST81 /N0An ðxÞS18 L0 In a similar manner 1 ½ 2 kn 88 1 ½ 3 kn 88 1 ½ 4 kn 88 2 ½ 1 kn 88 2 ½ 2 kn 88 2 ½ 3 kn 88 3 ½ 1 kn 88 3 ½ 2 kn 88 3 ½ 3 kn 88 3 ½ 4 kn 88 3 ½ 5 kn 88 4 ½ 1 kn 88 4 ½ 2 kn 88 4 ½ 3 kn 88 4 ½ 4 kn 88

References

¼ ½ 12 Hn /NAn ðxÞST81 /NBn ðxÞS18 L0 ¼ ½ 13 H/NAn ðxÞST81 /NCn ðxÞS18 L0 ¼ ½ 14 H/NAn ðxÞST81 /N00Cn ðxÞS18 L0 ¼ ½ 21 Hn /NBn ðxÞST81 /NAn ðxÞS18 L0 ¼ ½ 22 H/NBn ðxÞST81 /N0Bn ðxÞS18 L0 ¼ ½ 23 Hn /NBn ðxÞST81 /N0Cn ðxÞS18 L0 ¼ ½ 31 Hn /NCn ðxÞST81 /NAn ðxÞS18 L0 ¼ ½ 32 H/NCn ðxÞST81 /N00An ðxÞS18 L0 ¼ ½ 33 Hn /NCn ðxÞST81 /N0Bn ðxÞS18 L0 ¼ ½ 34 Hn /NCn ðxÞST81 /N0Cn ðxÞS18 L0 L ¼ ½ 35 H/NCn ðxÞST81 /N000 Cn ðxÞS18 0

¼ ½ 41 H/N0Cn ðxÞST81 /N0An ðxÞS18 L0 ¼ ½ 42 Hn /N0Cn ðxÞST81 /NBn ðxÞS18 L0 ¼ ½ 43 Hn /N0Cn ðxÞST81 /NCn ðxÞS18 L0 ¼ ½ 44 H/N0Cn ðxÞST81 /N00Cn ðxÞS18 L0

ðE:1a2pÞ

The total stiffness matrix contribution to mode n is then

½Kn 88 ¼

¼4 ¼ 3 jX 2pEhR iX i ½k  2 1n i ¼ 1 j ¼ 1 j n 88

2pEhR n 3 o 2pEhR n 1 o n 2 o n 3 o n 3 o f þ 6fn þ 6fn þ 8fn þ 9fn 1n2 6 n 1n2 n o 3 ðE:5Þ þ 10 f n

ðE:2Þ

[1] H. Ohtsubo, O. Watanabe, Stress analysis of pipe bends by ring elements, J. Pressure Vessel Technol. 100 (1978) 112–122. [2] H. Takeda, S. Asai, K. Kwata, A new finite element for structural analysis of piping systems, in: Proceedings of the Fifth SMIRT Conference, Berlin, 1979. [3] I.F. Pinto Correia, J.I. Barbosa, C.M. Mota Soares, C.A. Mota Soares, A finite element semi-analytical model for laminated axisymmetric shells: static, dynamic and buckling, Comput. Struct. 79 (2000) 299–317. [4] K.J. Bathe, C.A. Almeida, A simple and effective pipe elbow element—linear analysis, Trans. ASME J. Appl. Mech. 47 (1980) 93–100. [5] H. Santos, C.M. Mota Soares, C.A. Mota Soares, J.N. Reddy, A semi-analytical finite element model for the analysis of laminated 3D axisymmetric shells: bending, free vibration and buckling, Compos. Struct. 71 (2005) 273–281. [6] C. Militello, A.E. Huespe, A displacement-based pipe elbow element, Comput. Struct. 29 (1988) 339–343. [7] P.B. Goncalves, N.R.S.S. Ramos, Free vibration analysis of cylindrical tanks partially filled with liquid, J. Sound Vib. 195 (1996) 429–444. [8] E.M.M. Fonseca, F.J.M.Q. de Melo, C.A.M. Oliveira, Numerical analysis of piping elow for in-plane bending and internal pressure, Thin Walled Struct. 44 (2006) 393–398. [9] D. Zhou, Y.K. Cheung, S.H. Lo, F.T.K. Au, 3D vibration analysis of solid and hollow circular cylinders via Chebyshev–Ritz method, Comput. Methods Appl. Mech. Eng. 192 (2003) 1575–1589. [10] K. Weicker, R. Salahifar, M. Mohareb, Shell analysis of thin-walled pipes. Part II—finite element formulation, Int. J. Pressure Vessels Piping 87 (2010) 414–423. [11] S. Timoshenko, S. Woinowsky-Krieger, Theory of Plates and Shells, second ed., McGraw-Hill, New York, 1959. [12] E.M.M. Fonseca, F.J.M.Q. de Melo, C.A.M. Oliveira, Determination of flexibility factors in curved pipes with end restraints using a semi-analytical formulation, Int. J. Pressure Vessels Piping 79 (2002) 829–840. [13] S.A. Karamanos, Bending instabilities of elastic tubes, Int. J. Solid Struct. 39 (2002) 2059–2085.

92

R. Salahifar, M. Mohareb / Finite Elements in Analysis and Design 52 (2012) 83–92

[14] S.A. Karamanos, E. Giakoumatos, A.M. Gresnigt, Nonlinear response and failure of steel elbows under in-plane bending and pressure, J. Pressure Vessel Technol. 125 (2003) 393–402. [15] L. Xuebin, A new approach for free vibration analysis of thin circular cylindrical shell, J. Sound Vib. 296 (2006) 91–98. [16] L. Gan, X. Li, Z. Zhang, Free vibration analysis of ring-stiffened cylindrical shells using wave propagation approach, J. Sound Vib. 326 (2009) 633–646. [17] K. Weicker, R. Salahifar, M. Mohareb, Shell analysis of thin-walled pipes. Part I—field equations and solution, Int. J. Pressure Vessel Piping 87 (2010) 402–413. [18] R. Salahifar, M. Mohareb, Analysis of circular cylindrical shells under harmonic forces, Thin Walled Struct. 48 (2010) 528–539. [19] J. Mackerle, Finite elements in the analysis of pressure vessels and piping—a bibliography (1976–1996), Int. J. Pressure Vessels Piping 69 (1996) 279–339. [20] J. Mackerle, Finite elements in the analysis of pressure vessels and piping, an addendum (1996–1998), Int. J. Pressure Vessels Piping 76 (1999) 461–485. [21] J. Mackerle, Finite elements in the analysis of pressure vessels and piping, an addendum: a bibliography (1998–2001), Int. J. Pressure Vessels Piping 79 (2002) 1–26. [22] J. Mackerle, Finite elements in the analysis of pressure vessels and piping, an addendum: a bibliography (2001–2004), Int. J. Pressure Vessels Piping 82 (2005) 571–592. [23] K.J. Bathe, C. Almeida, A simple and effective pipe elbow element— interaction effects, J. Appl. Mech. 49 (1982) 165–171. [24] K.J. Bathe, C. Almeida, A simple and effective pipe elbow element—some nonlinear capabilities, Comput. Struct. 17 (1983) 659–667.

[25] M. Abo-Elkhier, Analysis of pipe bends using pipe elbow element, Comput. Struct. 37 (1990) 9–15. [26] F.J.M.Q. Melo, P.M.S.T. de Castro, The linear elastic stress analysis of curved pipes under generalized loads using a reduced integration finite element, J. Strain Anal. Eng. Des. 32 (1997) 47–59. [27] K.R. Sivadas, N. Ganesan, Dynamic analysis of circular cylindrical shells with material damping, J. Sound Vib. 166 (1993) 103–116. [28] A.J. Stanley, N. Ganesan, Dynamic response of cantilever cylindrical shells with discontinuity in thickness subjected to axisymmetric load, Comput. Struct. 55 (1995) 667–685. [29] Z.C. Xi, L.H. Yam, T.P. Leung, Semi-analytical study of free vibration of composite shells of revolution based on the Reissner–Mindlin assumption, Int. J. Solid Struct. 33 (1996) 851–863. [30] A.M. Yan, R.J. Jospin, D.H. Nguyen, An enhanced pipe elbow element—application in plastic limit analysis of pipe structures, Int. J. Numer. Methods Eng. 46 (1999) 409–431. [31] R. Salahifar, Analysis of Pipeline Systems under Harmonic Forces, Ph.D. Thesis, Department of Civil Engineering, University of Ottawa, Ottawa, ON, Canada, 2011. [32] W. Flugge, Stresses in Shells, Springer-Verlag, New York, 1973. [33] F.I. Niordson, Shell Theory, North-Holland Series In Applied Mathematics and Mechanics (1985). [34] V.V. Novozhilov, Thin Shell Theory, P. Noordhoff Ltd, Groningen—The Netherlands, 1964. [35] A.S. Saada, Elasticity: Theory and Applications, Pergamon Press Inc., New York, 1974.