Document Type : Research Paper
Authors
Composite Research Centre, Malek Ashtar University of Technology, Lavizan, Theran, Iran
Abstract
Keywords

Mechanics of Advanced Composite Structures 6 (2019) 201  223


Semnan University 
Mechanics of Advanced Composite Structures journal homepage: http://MACS.journals.semnan.ac.ir 
A New ThreeDimensional Refined HigherOrder Theory for Free Vibration Analysis of Composite Circular Cylindrical Shells
A. Davar^{*}, R. Kohandani, H.M. Panahiha
Composite Research Centre, Malek Ashtar University of Technology, Lavizan, Theran, Iran
Paper INFO 

ABSTRACT 
Paper history: Received 20180518 Received in revised form 20180827 Accepted 20190425 
A new closed form formulation of threedimensional (3D) refined higherorder shell theory (RHOST) to analyze the free vibration of composite circular cylindrical shells has been presented in this article. The shell is considered to be laminated with orthotropic layers and simply supported boundary conditions. The proposed theory is used to investigate the effects of the inplane and rotary inertias as well as transverse normal and shear strains on the dynamic response of thick composite cylindrical shells. The trapezoidal shape factor of the shell element is incorporated to obtain accurate stressresultants. Using Hamilton’s principle, the equations of motion are obtained and solved in terms of the Galerkin method. Numerical results for the natural frequencies are verified by making comparison with the 3D exact elasticity iterative solutions in the literature. In addition, the validity of the results is further verified by ABAQUS. According to the results, for thick composite cylinders with large lengthtoradius and orthotropic ratios, through thickness exact integration yields accurate stressresultants for proper prediction of the natural frequencies. 



Keywords: Free vibration Higherorder shear deformation theory Trapezoidal shape factor Circular cylindrical shells Composite 


© 2019 Published by Semnan University Press. All rights reserved. 
Cylindrical shells are widely used in many industries such as gas pipelines, petrol conveying. Also, cylindrical structures are common in modern industries such as aerospace, aircraft and marine structures. Based on classical shells theories, which are based on KirchhoffLove’s hypothesis, many studies have been performed on shells [1]. Although classical shell theories ignore the transverse stress and strain components for easy calculation, this omission gives inadequate results for the analysis of thick cylindrical shells [1]. Some research studies are presented in the literature that investigate the effects of shear deformation for dynamic response of composite cylindrical shells [2]. Leissa [3] has summarized many studies in the stateoftheart in his research work. According to these research studies, the effect of shear deformation can become significant for small lengthtothickness or radiustothickness ratios. Bhimaraddi [4], developed a twodimensional higherorder shell theory to investigate the dynamic response of composite circular cylindrical shell and the traction free condition is assumed for inner and outer surfaces of the shell. Reddy and Liu [5] presented a twodimensional (2D) higherorder theory for laminated elastic shells. The theory accounts for parabolic distribution of the transverse shear strains through thickness of the shell and tangential stressfree boundary conditions on the boundary surface of the shell.
The 2D higherorder shell theories consider the effects of shear deformation and rotary inertia and they are more useful than the thin shell theories for the analysis of moderately thick shell structures. In order to analyze the thick shells, 2D higherorder shell theories are not adequate especially in the case of higher frequencies. In order to analyze the thick shells, the transverse normal stress and strain components which are neglected in the 2D higherorder shell theories, should be accounted for in the analysis which is based on threedimensional (3D) shell theories.
Due to accounting all the transverse stress and strain components (which are ignored in the 2D higherorder shell theories), the dynamic analysis of circular cylindrical shells on the basis of the governing equations of the 3D elasticity attracted the attention of researchers. In recent years, by refinement of thickshell theories, some new 3D shell theories for the case of homogeneous cylindrical shells were investigated [68] as reviewed by Qatu [9]. Khalili et al. [10] investigated dynamic responses of free vibration analysis of homogenous isotropic circular cylindrical shells based on a new 3D refined higherorder theory.
In the case of multilayered anisotropic composite shells, the effects of transverse shear deformation are more significant as compared to isotropic shells. Hence, the dynamic behavior of composite shells is more complicated than isotropic ones. Because of this complexity, accurate results for dynamic response of composite shells need threedimensional modeling instead of twodimensional one especially for analysis of thick shells where transverse normal and shear strains become more significant. Rogers and Knight [11] have formulated a linear higherorder finite element to analyze an axisymmetric composite structure. A higherorder theory for the analysis of composite cylindrical shells was proposed by Murthy et al. [12] by expanding the displacement variables in the form of power series and retaining a finite number of terms. As a result, the formulation allows for arbitrary variation of inplane displacement. Threedimensional elasticity solutions were presented for the vibration of crossply laminated simply supported cylindrical shells by Ye and Soldatos [13]. They used an iterative procedure and after a few iterations, they obtained the exact values for the natural frequencies. Natural frequencies and their mode shapes of some homogeneous orthotropic crossply cylinders were investigated. Kant and Menon [14] presented a higherorder refined theory for composite and sandwich cylindrical shells with finite number of elements which is suitable for the analysis of thin and moderately thick anisotropic laminated cylindrical shells. Timarci and Soldatos [15] presented comparative dynamic research studies for symmetric crossply cylindrical shells using unified sheardeformable shell theory.
Most of the research studies for higherorder shear deformation theories that include shear deformation and rotary inertia, failed to consider the ( terms (trapezoidal shape factor) that is considered due to the fact that the stresses over the thickness of the shell have to be integrated on a trapezoidal crosssection of a shell element to obtain the accurate stress resultants. As shown in Fig. 1, an element of the shell section is presented. As it can be seen, taking into account the large shape trapezoidal coefficient (including 1+z/R terms), instead of the rectangular shape (excluding 1+z/R terms), is closer to reality, and therefore precision of the integration increased for calculating the stress resultants in the axial direction.
Chang [16] and Leissa and Chang [17] considered this term but by neglecting the terms beyond the order of . For the first time, Qatu [18] utilized the ( shape factor within the framework of first order shear deformation theory (FSDT) to analyze the free vibration of laminated deep thick shells. Lam and Qian [19] developed a theoretical analysis and analytical solution for vibrations of thick symmetric angleply laminated composite shells considering trapezoidal shape factor ( . Icardi and Ruotolo [20] presented a multilayered model based on a secondorder expansion of the ( terms. They presented some numerical results concerning eigenfrequencies and stress distributions across the thickness of simply supported, crossply cylindrical shells. As a result, incorporation of the secondorder expansion of the ( terms appears to be suited for technical purposes, as it can improve the accuracy for predicting the overall and local behavior of rather thick shells. Other research studies which incorporate the ( terms in the static and dynamic analysis of thick laminated cylindrical shells are presented in Refs. [2125]. In these studies, the most popular procedures are finite element method, Ritz method and the series solution method.
Section Development 
(a) (b)
Fig. 1. (a) Circumferential crosssection of a thick cylinder for integrating stress resultants in axial direction; (b) Regions shown by + and  signs indicate the area differences between the assumed (rectangular and trapezoidal) developed shapes of the crosssection.
The main purpose of this work is to investigate a closed form solution of the free vibration of simply supportedsimply supported (SSSS) composite laminated circular cylindrical shells using a threedimensional (3D) refined higherorder shell theory (RHOST). The effects of the inplane and rotary inertias and transverse normal and shear strains on the dynamic response of composite cylindrical shells have been investigated. Due to the fact that the stresses over the thickness of the shell are to be integrated on trapezoidallike cross section of a shell element, trapezoidal shape factor ( is also considered for the first time in the framework of the present RHOST. The present work is an extension of the first author earlier research on free vibrations of thick homogenous isotropic cylinders [10] to multilayered thick composite cylinders. The advantage of the present RHOST is that no iterative procedure like those used for example in Refs. [1] and [13] is required for calculating the natural frequencies and hence, less CPUtime is consumed and this would be useful especially in optimization processes where frequency should be calculated several times.
A circular cylindrical shell as shown in Fig. 2 is considered with radius R, thickness h and length L. The displacement components in the axial, tangential and radial directions are u, v and w, respectively and the reference coordinate system (x,j, z), is placed on the middle surface of the cylindrical shell.
Fig. 2. A circular cylindrical shell with the reference coordinate system
In order to formulate a 3D elasticity problem, the Taylor’s series expansion is used and the following equations are obtained by expanding the displacement components u(x,j, z, t), v(x,j, z, t) and w(x,j, z, t) in terms of thickness coordinate z of any point of shell space [10]:
(1) 
The terms u, v and w are the displacements components and t is the time. It should be noted that 12 displacement parameters are presented in Eq. (1) as a higherorder displacement field. By setting the coefficient equal to 1 in Eq. (1), the trapezoidal shape factor of the cylindrical shell is applied in the equilibrium equations and the HOST12 theory ( is refined to RHOST12. , are the inplane displacements of the cylindrical shell and is the transverse displacement of a point (x,j) on the shell middle surface. , are the rotation functions of the normal to the shells middle surface about j and x axis, respectively. , , , , and are the higherorder terms in the Taylor’s series expansion that represent higherorder transverse deformation modes. For the firstorder shear deformation theory, only , , , and are considered as displacement filed. The general straindisplacement relations in the cylindrical coordinate system according to linear theory of elasticity for circular cylindrical shells are defined as follow [10]:
(2) 
By substituting Eq. (1), the expressions for displacement at any point within the shell, the linear strains in terms of middle surface displacements are obtained as follow:
(3) 
where:
(4)) 

For an orthotropic material, 3D stressstrain relations are obtained by Hooke’s law as [1]:
(5) 
coefficients are defined as:
(6) 
where are Young’s modulus of elasticity, are Poisson’s ratio, are the shear moduli for composite material in different directions. The relation between offaxis stress and strain for the layer of a multilayered composite cylindrical shell is defined as follows:
(7) 
where Q_{ij} are elements of the reduced stiffness matrix as defined in Appendix A.
By substituting Eq. (3) into Eq. (7) and integrating through the shell thickness, Eq. (7) is refined to:
(8) 

(9) 
The matrices , and are given in Appendix B. is the shear correction factor whose value is considered equal to 1 for higherorder theories and equal to first order shear deformation theory. In Appendix C, accurate method of calculation of integrals thorough the shell thickness in stressresultant equations, including the ( terms, are presented. and , the middle surface strain vector and stressresultant vector, respectively, in Eq. (8) are defined as follow:
(10) 

(11) 
The components of the stressresultant vector for the composite shell are defined as:
(12) 
where NL is the number of composite layers.
Using Hamilton’s principle, the equations of motion for the free vibration analysis are obtained. It could be defined as follows in analytical form:
(13) 
where U is the total strain energy due to deformation, W is the potential of the external loads and K is the kinetic energy. Due to the assumption of the absence of damping and external loads, the Hamilton’s principle could be summarized as follows:
(14) 
U, the total strain energy due to deformationin Eq. (13) is defined as:
(15) 
and the kinetic energy, K is defined as:
(16) 
ρ is the mass density of the material of the shell and ( represents differentiation with respect to time. In Eqs. (15) and (16), the differential element on the middle surface of the shell is defined as [10]:
(17) 
Substituting the appropriate strain expressions given by Eq. (4) and the displacement expressions given by Eq. (1) in Eq. (14) and integrating the resulting expressions by parts, after separating the coefficients of , , , , , , , , , , , , the equations of motions are obtained:

(18a) 


(18b) 

(18c) 

(18d) 

(18e) 

(18f) 

(18g) 

(18h) 

(18i) 

(18j) 

(18k) 

(18l) 

where the inertia terms in the right side of Eqs. (18) are given by
(19) 
Natural and essential boundary conditions for simply supported conditions at x=0 and x=L are as follows:
(20) 
The Galerkin method is used to solve the free vibration problem of simply supportedsimply supported circular cylindrical composite shells. The boundary conditions for simply supported edges at x=0 and x=L is applied as follows [10]:
(21) 
In order to satisfy the boundary conditions, the displacement components are expanded as follow [6]:
(22) 
where in Eq. (22), = and are the natural angular frequencies (in ) and are related to the mode numbers (m,n) where m is the axial halfwave number and n is the circumferential wave number. , , , , , , , , , , , are the constant amplitudes of vibrations related to the natural mode shapes. By substituting Eq. (22) into Eqs. (18) and applying the Galerkin method, after simplification and collecting coefficients, the following eigenvalue equation is obtained:
(23) 
where and d is the displacement vector, corresponding to the mode shape numbers (m,n). Generally, between the 12 eigenvalues (frequencies) obtained from Eq. (23), the smallest one is associated to the bending vibration mode shape corresponding to the specified mode numbers (m,n). The lowest eigenvalue is called fundamental frequency of bending vibration. The elements of the stiffness matrix [K] and the mass matrix [M] are given in Appendix D.
In order to analyze the free vibration of composite circular cylindrical shells with simply supportedsimply supported boundary conditions, a computer code using MATLABR13 based on the formulation of the present shell theories is developed. Different examples of composite cylindrical shells with a wide range of thicknesstoradius ( and lengthtoradius ( ratios are investigated to show the efficiency and accuracy of the present formulations. In order to verify the present results, they were compared to the analytical results available in the literature. Furthermore, the results were validated with those obtained using Lanczos eigenfrequency extraction subroutine in ABAQUS/Standard code. In order to obtain accurate results of the free vibration analysis, the stress resultants were calculated using exact integration over the thickness of the composite cylindrical shells. In addition, in contrast to some 3D elasticity theories, e.g. Refs. [1, 13], in the literature for the free vibration analysis, the present RHOST does not require any iterative procedure and convergence study. This is an advantage from the sense that computational time in the present RHOST is less than these iterative procedures.
Unless otherwise stated, the following geometrical and material properties are used hereinafter:
(24) 
Also, unless otherwise mentioned, the natural frequency parameter is considered to be:
(25) 
For the verification purpose, in Tables 1 to 4, the results of the present RHOST12 and HOST12 are compared to those obtained from 3D elasticity exact solution method, reported by Ref. [13].
In Table 1, the lowest natural frequency parameters, for 4layered cylindrical shells having symmetric crossply layup are presented. As can be seen in this Table, in all cases, the first frequency parameters of layup are greater than those for layup. Furthermore, good accuracy is obtained in comparison with the results of Ref. [13] for different values of thickness ratios and mode shape numbers. By increasing h/R, the discrepancies are increased. Also, by increasing n, the discrepancies are increased for [0/90]_{s} layup, unless at n=3 for . The maximum discrepancy (1.75%) is corresponded to h/R=0.3 and n=2 for [90/0]_{s} layup.
In Tables 2 and 3, the first three frequency parameters of crossply composite cylindrical shells are shown for different values of , and circumferential mode number (n). The results were also validated by making comparison with Ref. [13]. According to these Tables, the natural frequency parameters of the laminated cylinders increased by increasing the number of layers. This trend occurs due to the fact that by increasing the number of layers, the bendingextensional coupling decreased for an antisymmetric crossply laminate. It is necessary to be noted that in Tables 2 and 3, by increasing both h/R and n, the discrepancies increased. The maximum discrepancy for the present RHOST12 (2.84%) and for the present HOST12 (5.12%) are both corresponded to h/R=0.3 and n=3 for [0/90] layup in Table 2. As shown in these Tables, all the discrepancies corresponding to the present RHOST12 are very less than those for the present HOST12. This outcome reveals the importance of incorporating the trapezoidal shape factor (1+z/R terms) in the present analytical formulations. In addition, as can be seen in these Tables, except for h/R=0.1 for layup [0/90]_{2} in Table 2, generally by increasing the number of layers, the discrepancies decreased. Also, good agreement between the present results of RHOST12 and Ref. [13] results shows the accuracy of the present theory.
Table 4 shows the values of the first three frequency parameters for different orthotropic ratios of 2layered unsymmetric crossply composite cylindrical shells for different values of the thicknesstoradius ratio . The present RHOST12 results were validated by making comparison with those reported by Ref. [13] and good agreement was observed. According to the results, all frequency parameters increased by increasing either the stiffness ratio or the thicknesstoradius ratio of the cylinders. Also, the discrepancy increased by increasing the values of and h/R. The maximum discrepancy (5.6%) is corresponded to =40 and h/R=0.5 for the second (II) frequency. There is not a specific trend for the discrepancies when the frequency number (I, II or III) is changed.
In Fig. 3, variations of the natural frequency parameter vs. thicknesstoradius ( is indicated for symmetric crossply layups. Also, the results of the present RHOST12 and HOST12 theories for L/R= 10 are compared with the present FEM analysis. According to Fig. 3(a) by decreasing the value of orthotropic ratio, , the discrepancies between the results of the present RHOST12 and HOST12 and the present FEM results increased specially for greater values of .
Table 1. Natural frequency parameters, , for composite circular cylindrical shells with symmetric crossply layups (m=1)
h/R 
Theory 
[0/90]_{s} 

[90/0]_{s} 



n=1 
n=2 
n=3 

n=1 
n=2 
n=3 
0.1 
RHOST12 (present) 
0.079302 
0.066377 
0.064700 

0.070809 
0.052872 
0.059267 
Ref. [13] 
0.079277 
0.066335 
0.064600 

0.070738 
0.052748 
0.059130 


0.03^{*} 
0.06 
0.15 

0.10 
0.23 
0.23 










0.2 
RHOST12 (present) 
0.175333 
0.163123 
0.171778 

0.151538 
0.131548 
0.160240 
Ref. [13] 
0.175188 
0.162844 
0.170868 

0.150651 
0.130168 
0.158886 


0.08 
0.17 
0.53 

0.58 
1.06 
0.85 










0.3

RHOST12 (present) 
0.273215 
0.263860 
0.286369 

0.239027 
0.222623 
0.272072 
Ref. [13] 
0.272860 
0.263048 
0.283798 

0.236385 
0.218779 
0.268258 


0.13 
0.30 
0.90 

1.11 
1.75 
1.42 
^{ *}Percentage discrepancy ((Present –Ref. [13])/ Ref. [13])*100.
Table 2. Natural frequency parameters, , for composite circular cylindrical shells with unsymmetric crossply layups (m=1)
h/R 
Theory 














n=1 

n=2 

n=3 


n=1 

n=2 

n=3 

0.1 
RHOST12 (present) 
0.069519 
0.13^{*} 
0.049802 
0.35 
0.046207 
0.56 

0.074085 
0.22 
0.058276 
0.52 
0.059502 
0.78 
HOST12 (present) 
0.069594 
0.24 
0.049986 
0.72 
0.046627 
1.46 

0.074141 
0.30 
0.058461 
0.84 
0.059881 
1.42 

Ref. [13] 
0.069428 

0.049630 

0.045949 


0.073919 

0.057975 

0.059043 

















0.2 
RHOST12 (present) 
0.148051 
0.84 
0.122233 
1.64 
0.130875 
1.99 

0.162243 
0.81 
0.145715 
1.48 
0.162551 
1.77 
HOST12 (present) 
0.148539 
1.17 
0.123520 
2.72 
0.133160 
3.77 

0.162664 
1.08 
0.146790 
2.23 
0.164226 
2.82 

Ref. [13] 
0.146819 

0.120255 

0.128317 


0.160932 

0.143589 

0.159729 

















0.3 
RHOST12 (present) 
0.233711 
1.61 
0.208223 
2.64 
0.232943 
2.84 

0.253687 
1.10 
0.239751 
1.82 
0.272722 
2.01 
HOST12 (present) 
0.235187 
2.24 
0.211727 
4.37 
0.238106 
5.12 

0.254902 
1.59 
0.242367 
2.94 
0.276154 
3.29 

Ref. [13] 
0.230019 

0.202861 

0.226517 


0.250922 

0.235457 

0.267347 

^{*}Percentage discrepancy ((Present –Ref. [13])/ Ref. [13])*100.
Table 3. Natural frequency parameters, for composite circular cylindrical shells with unsymmetric crossply layups (m=1)
h/R 
Theory 














n=1 

n=2 

n=3 


n=1 

n=2 

n=3 

0.1 
RHOST12 (present) 
0.075019 
0.11^{*} 
0.059930 
0.24 
0.062060 
0.36 

0.075390 
0.07 
0.060569 
0.15 
0.063045 
0.24 
HOST12 (present) 
0.075052 
0.15 
0.060071 
0.48 
0.062379 
0.79 

0.075408 
0.09 
0.060679 
0.33 
0.0632470 
0.56 

Ref. [13] 
0.074939 

0.059787 

0.061838 


0.075339 

0.060477 

0.062896 

















0.2 
RHOST12 (present) 
0.165499 
0.39 
0.151194 
0.72 
0.170352 
0.90 

0.166894 
0.27 
0.153521 
0.49 
0.173719 
0.64 
HOST12 (present) 
0.165812 
0.58 
0.152001 
1.26 
0.171439 
1.55 

0.167122 
0.41 
0.154138 
0.89 
0.174432 
1.05 

Ref. [13] 
0.164852 

0.150114 

0.168829 


0.166445 

0.152779 

0.172616 

















0.3 
RHOST12 (present) 
0.259425 
0.57 
0.249229 
0.99 
0.285798 
1.20 

0.262183 
0.42 
0.253746 
0.75 
0.292216 
0.99 
HOST12 (present) 
0.260390 
0.95 
0.251165 
1.78 
0.288778 
1.93 

0.262929 
0.71 
0.255779 
1.32 
0.293348 
1.38 

Ref. [13] 
0.257947 

0.246783 

0.282406 


0.261081 

0.251849 

0.289353 

^{*}Percentage discrepancy ((Present –Ref. [13])/ Ref. [13])*100.
Table 4. First three lowest frequency parameters, for composite circular cylindrical shells with unsymmetric crossply layups (m=n=1)

Theory 


h/R=0.1 

h/R=0.3 

h/R=0.5 



I 
II 
II 

I 
II 
II 

I 
II 
II 
10 
RHOST12 (present) 
0.06192 
0.15824 
0.29462 

0.20945 
0.47465 
0.71821 

0.38460 
0.78261 
0.98567 
Ref. [13] 
0.06192 
0.15824 
0.29444 

0.20878 
0.47432 
0.71297 

0.38198 
0.78053 
0.96867 


0^{*} 
0 
0.06 

0.32 
0.06 
0.73 

0.68 
0.26 
1.75 














20 
RHOST12 (present) 
0.06632 
0.20305 
0.38953 

0.22237 
0.58876 
0.80612 

0.39948 
0.92867 
1.06172 
Ref. [13] 
0.06629 
0.20302 
0.38888 

0.22057 
0.58695 
0.79149 

0.39425 
0.91408 
1.02944 


0.04 
0.01 
0.16 

0.81 
0.3 
1.84 

1.32 
1.59 
3.13 














30 
RHOST12 (present) 
0.06828 
0.23929 
0.45417 

0.22921 
0.66924 
0.85234 

0.40660 
1.00773 
1.10908 
Ref. [13] 
0.06823 
0.23922 
0.45285 

0.22638 
0.66463 
0.82864 

0.39975 
0.97244 
1.07140 


0.07 
0.02 
0.29 

1.25 
0.69 
2.86 

1.71 
3.62 
3.51 














40 
RHOST12 (present) 
0.06951 
0.27021 
0.50185 

0.23371 
0.72855 
0.88451 

0.41091 
1.05494 
1.14304 
Ref. [13] 
0.06943 
0.27009 
0.49971 

0.23002 
0.71976 
0.85292 

0.40298 
0.99892 
1.10068 


0.11 
0.04 
0.42 

1.6 
1.22 
3.7 

1.96 
5.6 
3.84 
^{*}Percentage discrepancy ((Present –Ref. [13])/ Ref. [13])*100.
(a) 
(b) 
Fig. 3. Effect of layup on the lowest natural frequency parameter, , SSSS composite cylinder vs. thicknesstoradius ratio ( ); . (a). =40 (L/R=1). (b). =1 (L/R=10) 

As illustrated in Fig. 3(a), at h/R=1.9, the discrepancies of the present RHOST12 and HOST12 for =40 are 0.23% and 6.21% for layup and 3.19% and 6.66% for layup, respectively. At =1.9, the discrepancies of the present RHOST12 and HOST12 for =1 are 0.44% and 18.4% for layup and 0.43% and 19.33% for layup, respectively, as indicated in Fig. 3(b). Generally, the frequencies corresponded to layup are greater than layup.
In Figs. 4(a) and 4(b), variations of the lowest natural frequency parameter vs. ratio is indicated for [0/90]s and [0/90]2 layups of a crossply composite circular cylindrical shell for the present RHOST12 and HOST12. According to these figures, regardless of the layup sequence (symmetric or unsymmetric) for and 1.5, by increasing ratio from 5 to 20, the differences between the present RHOST12 and HOST12 increased from about 0% and 2% to about 2.63% and 6.7%, respectively.
Variations of the lowest natural frequency parameter, vs. ratio, for different orthotropic ratios ( are presented in Fig. 5. The results of the present HOST12 and RHOST12 are compared to each other for =1 (Fig. 5(a)) and also to the present FEM results for =10 (Fig. 5(b)).
(a) 
(b) 
Fig. 4. Lowest natural frequency parameter, , SSSS composite cylinder vs. lengthtoradius ratio ( ) for two different layups ; (a). =0.5 ( =1). (b). h/R=1.5 ( =1). 

(a) 
(b) 
Fig. 5. Lowest natural frequency parameter for SSSS composite cylinder vs. thicknesstoradius ratio ( ) for different orthotropic ratios ( ; (a). =1, (b). L/R=10. 
As it can be seen in Fig. 5(a), regardless of the value of , the difference between the present RHOST12 and HOST12 results are not considerable since the shell length is short. However, according to Fig. 5. (b), by increasing from 0.1 to 1.9, the discrepancies between the present RHOST12 and FEM for =1, 10 and 40 increased from about 0% to about 0.67% 3.35% and 1.64%, and these discrepancies for HOST12 increased from about 0% to about 18.58% , 17.87% and 15.7%, respectively.
According to Fig. 5(b), there is a good agreement between the present RHOST12 and FEM results and noticeable discrepancies were found between HOST12 and FEM results.
Fig. 6(a) indicates the effect of different layups on the lowest natural frequency parameter vs. h/R ratio for L/R=1. In this figure, the maximum difference between the results of the present RHOST12 and HOST12 is about 1.01%. However, in the case of L/R=10, as indicated in Fig. 6(b), the geometric parameter L/R has considerable influence on the accuracy of the present theories. By increasing the value of L/R, the differences between the present HOST12 and RHOST12 increased. The results of the present FEM simulations are also compared. As it can be seen in Fig. 6(b), the maximum discrepancy 8.5% occurs in the case of [90] layup between the present HOST12 and FEM results.
It could be observed from Figs. 6(a) and 6(b) that for both RHOST12 and HOST12 by increasing the volume fraction of zero angle layers in the laminate, the frequency of the cylinder increased.
Figs. 7(a) and 7(b) illustrate the variations of frequency parameter vs. orthotropic ratio E1/E2 for L/R=1 and 10, respectively. In Fig. 7(a), a small difference between the present RHOST12 and HOST12 exists and it is almost unchanged by increasing the value of E1/E2 the present RHOST12 and HOST12. In addition, the frequency converges to an almost constant value. However, in case of L/R=10 in Fig. 7(b), by increasing E1/E2 and h/R, the differences between the present RHOST12 and HOST12 increased.
The maximum discrepancy (14.55%) between the present theories and the present FEM simulations occurs for HOST12 at h/R=1.8 and E1/E2 =45 as it can be seen in Fig. 7(b). In fact, by increasing the value of L/R, the influence of the exact integration of the stress resultants over the trapezoidallike crosssection of the shell becomes more important. Since in the present HOST12, this important point is not considered, this theory fails to predict the correct values of the frequency in contrast to the present RHOST12 especially for higher values of h/R, E1/E2 and L/R ratios.
(a) 
(b) 
Fig. 6. Lowest natural frequency parameter for SSSS composite cylinder vs. thicknesstoradius ratio ( ) for different crossply layups. (a). =1. (b). =10 
(a) 
(b) 
Fig. 7. Lowest natural frequency parameter for SSSS composite cylinder vs. orthotropic ratio ( ) for [0/90] layup. (a). =1. (b). L/R= 10. 
In Fig. 8, variations of lowest natural frequency parameters, ω^*, vs. number of layers (n) in [0/90]n layup for different orthotropic ratios (E_1⁄E_2 ) have been investigated. As shown in Fig. 8(a), for L/R=1 and h/R=0.1, regardless of the value of E1/E2, by increasing the number of layers (n) in [0/90]n layup, no considerable difference could be observed between the present RHOST12 and HOST12 results. Also, by increasing the number of layers in [0/90]n layup, regardless of the value of E1/E2, the frequency converges to a special constant value. Furthermore, in Fig. 8(b), for L/R=10 and h/R=1.5, the frequency converges to another special constant value. However, in contrast to Fig. 8(a), considerable difference could be observed between the present RHOST12 and HOST12 results by changing the value of E1/E2. For E1/E2=1, the differences between the present RHOST12 and HOST12 are almost unchanged by increasing the number of layers (n) in [0/90]n layup. While, for E1/E2=10 and 40, there is a special value for the number of layers (n) where the difference between the present RHOST12 and HOST12 is negligible and before and after this special value, this difference becomes clear especially for lower values of the number of layers (n).
In Table 5, the lowest natural frequency parameters, obtained from the present analytical theories are compared with those reported by Ref. [13]. In addition, results are compared to those obtained using Lanczos eigenvalue extraction method in ABAQUS/Standard solver and the associated mode shapes are depicted in Table 5. For the finite element 3D (FE) modeling of the composite cylindrical shells, 8noded continuum shell (SC8R) was used and convergence study for the elements size was achieved.
As it can be seen from Table 5, for the first bending mode No. (1,1), the absolute values of the discrepancies between the present theories and those of Ref. [13] for different thicknesstoradius ratios = 0.1, 0.2 and 0.3 are 0.13%, 0.64% and 1.61%, respectively, for RHOST12 and 0.24%, 1.17%and 2.24%, respectively, for HOST12. According to Table 5, for the second bending mode No. (1,2), the absolute values of discrepancies between the present theories and those of Ref. [13] for different thicknesstoradius ratios =0.1, 0.2and 0.3 are 0.35%, 1.64% and 2.64%, respectively, for RHOST12 and 0.72%, 3.96% and 7.44%, respectively, for HOST12.
Also, as shown in Table 5, for the third bending mode No. (1,3), the absolute values of the discrepancies between the present theories and those of Ref. [13] for different thicknesstoradius ratios = 0.1, 0.2and 0.3 are 0.56%, 1.99% and 2.84%, respectively, for RHOST12 and 1.41%, 5.43% and 8.87%, respectively, for HOST12. Hence, by increasing the mode number, generally the discrepancies increased for both the present RHOST12 and HOST12.
(a) 
(b) 
Fig. 8. Lowest natural frequency parameter for SSSS composite cylinder vs. number of layers (n) in [0/90]n layup for different orthotropic ratios ( ). (a). L/R=1 and h/R=0.1. (b). L/R=10 and h/R=1.5 
According to Table 5, the present FEM analysis also indicates good accuracy as compared to those results in Ref. [13]. For different thicknesstoradius ratios =0.1, 0.2and 0.3, for the first mode (1,1), the discrepancies between the present FEM results and Ref. [13] are 0.04%, 1.09% and 0.69%, respectively. For the second bending mode (1,2), the discrepancies are 0.02%, 0.11% and 1.43%, respectively, and for the third bending mode (1,3), the discrepancies are 0.24%, 1.78% and 3.47%, respectively. As compared to Ref. [13], in most cases, the discrepancies of the present FEM results are less than those for the present RHOST12. However, in some cases like mode No. (1,1) with = 0.2 and mode no. (1,3) with = 0.3, the results of the present RHOST are closer than the present FEM results to those for Ref. [13].
Table 5. Comparison of lowest natural frequency parameters for composite circular cylindrical shells with unsymmetric crossply layups ([0/90]) with those obtained using Lanczos method of eigenvalue extraction in ABAQUS/Standard solver and associated mode shapes (L/R=1).
Mode no. (m,n) 
Theory 
=0.1 

=0.2 

=0.3 

(1,1) 
Ref. [13] 
0.069428 

0.146819 

0.230019 

RHOST12 (present) 
0.069519 
0.13^{*} 
0.148051 
0.84 
0.233711 
1.61 

HOST12 (present) 
0.069594 
0.24 
0.148539 
1.17 
0.235178 
2.24 

FEM (present) 
0.069458 
0.04 
0.148423 
1.09 
0.231612 
0.69 

FSDT(present) 
0.069858 
0.61 
0.150135 
2.25 
0.239457 
4.1 

Mode shape 

(1,2) 
Ref. [13] 
0.049630 

0.120255 

0.202861 

RHOST12 (present) 
0.049802 
0.35 
0.122233 
1.64 
0.208223 
2.64 

HOST12 (present) 
0.049986 
0.72 
0.123520 
2.72 
0.211727 
4.37 

FEM (present) 
0.049618 
0.02 
0.120119 
0.11 
0.199962 
1.43 

FSDT(present) 
0.050340 
1.43 
0.125735 
4.55 
0.217378 
7.15 

Mode shape 

(1,3) 
Ref. [13] 
0.045949 

0.128317 

0.226517 

RHOST12 (present) 
0.046207 
0.56 
0.130875 
1.99 
0.232943 
2.84 

HOST12 (present) 
0.046620 
1.46 
0.133160 
3.77 
0.238106 
5.12 

FEM (present) 
0.045838 
0.24 
0.126032 
1.78 
0.218653 
3.47 

FSDT(present) 
0.047055 
2.40 
0.135579 
5.65 
0.244040 
7.73 

Mode shape 
^{*}Percentage discrepancy ((Present –Ref. [13])/ Ref. [13])*100.
For the first time, a closed form solution method for free vibration analysis of composite thin and thick simply supported cylindrical shells on the basis of 3D refined higherorder shell theory (RHOST) is presented in this study. The effect of the trapezoidal shape factor (1+z/R terms) of the crosssection of the orthotropic composite circular cylindrical shells wa incorporated exactly in the formulations. The characteristic eigenvalue equation was obtained based on Hamilton’s principle and by applying Galerkin method to the governing equations, natural frequencies were obtained. The applicability and validity of the present theory were confirmed by verifying the results with those obtained using the exact 3D elasticity method for a wide range of thicknesstoradius and thicknesstolength ratios. Comparisons of the results for thick cylindrical shells with published results in the literature were carried out and good agreement was observed.
The present theory does not require any convergence study, in contrast to some existing iterative approaches in the literature that require a few iterations to achieve sufficient convergence to the exact solution. This is an important advantage of the present RHOST. Furthermore, the natural frequencies associated to highermodes of moderately thick, thick and very thick composite cylinders, never published in the literature before, were compared to those obtained using FE modeling in ABAQUS commercial software. The results show that considering the effect of the (1+z/R) terms in the calculation of stress resultants, would lead to a reliable higherorder theory for the free vibration analysis of highly orthotropic composite circular cylindrical shells, especially for the cases of long and thick hollow cylinders.
Nomenclature
Shell’s thickness 
h 
Length of the cylinder 
L 
Mean radius of the cylindrical shell 
R 
Position coordinate in axial direction 
x 
Position coordinate in radial direction 
z 
Position coordinate in tangential direction 

Displacement component in axial direction 
u 
Displacement component in tangential direction 
v 
Displacement component in radial direction 
w 
Coefficient of trapezoidal shape 

Displacement component 

Strain components 

Vector of strain components 

Vector of stress resultants components 

Normal and shear stresses of each layer 

Normal and shear strains of each layer 

Elements of stiffness matrix 

Elements of reduced stiffness matrix 

Young's modulus 

Shear modulus 

Poisson coefficients 

The rotation angle of the fiber relative to the main axis 

Normal and shear stresses for a multilayer 

Number of layers 
NL 
Shell stiffness matrices 

Total strain energy 
U 
Energy from external forces 
W 
Total kinetic energy 
K 
Shell mass inertia 

Differential operators 

Number of halfaxial waves 
m 
Number of circumferential waves 
n 
Time functions in generalized coordinates 

Natural frequency for mode (m, n) 

Fundumental natural frequency 

Stiffness matrix 
K 
Constant natural mode shapes 

Mass matrix 
M 
Shear Correction Coefficient (first order shear deformation theory) 

Stress resultants 

Strain and curvature components 
Appendixes
Appendix A. Elements of Reduced Stiffness Matrix Q_{ij}
where and ; is fibre orientation (in radians) with respect to xaxis of the shell.
The terms , and (j=1, 2, … ,7), used in the following matrices ( , , , , and ), are defined in Appendix C.
=
In matrices , , , , and the terms , and are defined as follows:
C.1. Definition of
(C‑1) 
, 
C.2. Definition of
(C‑2) 

where are defined in Eq. (C1)
C.3. Definition of
(C‑3) 
In the case of :
(C‑4) 
In the case of after taking exact integration through the thickness of each layer, the following results were obtained:
(C‑5) 


D.1. Elements of Stiffness Matrix :
, , ,
, , , , , , , ,
, ,
, , ,
, , ,
, , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,
, , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,
, , , , , , , , , , ,
D.2. Elements of Mass Matrix :
, , , ; , , , ; , , , ; , , , ; , , , ; , , , ; , , , ; , , , ; , , , ; , , , ; , , , ; , , , ;
Other elements of M are equal to zero.
[1] Loy C, Lam K. Vibration of thick cylindrical shells on the basis of threedimensional theory of elasticity. Journal of sound and Vibration 1999; 226(4): 71937.
[2] Hildebrand F, Reissner E, Thomas G. Notes on the foundations of the theory of small displacements of orthotropic shells: National Advisory Committee for Aeronautics Washington, DC; 1949.
[3] Leissa AW. Vibration of shells: Scientific and Technical Information Office, National Aeronautics and Space Administration Washington, DC, USA; 1973.
[4] Bhimaraddi A. A higher order theory for free vibration analysis of circular cylindrical shells. International Journal of Solids and Structures 1984; 20(7): 62330.
[5] Reddy J, Liu C. A higherorder shear deformation theory of laminated elastic shells. International Journal of Engineering Science 1985; 23(3): 31930.
[6] Singal R, Williams K. A theoretical and experimental study of vibrations of thick circular cylindrical shells and rings. Journal of Vibration and Acoustics 1988; 110(4): 5337.
[7] Soldatos K, Hadjigeorgiou V. Threedimensional solution of the free vibration problem of homogeneous isotropic cylindrical shells and panels. Journal of Sound and Vibration 1990; 137(3): 36984.
[8] Kant T, Kumar S, Singh U. Shell dynamics with threedimensional degenerate finite elements. Computers & structures 1994; 50(1): 13546.
[9] Qatu MS. Recent research advances in the dynamic behavior of shells: 1989–2000, Part 2: Homogeneous shells. Applied Mechanics Reviews 2002; 55(5): 41534.
[10] Khalili S, Davar A, Malekzadeh Fard K. Free vibration analysis of homogeneous isotropic circular cylindrical shells based on a new threedimensional refined higherorder theory. International Journal of Mechanical Sciences 2012; 56(1): 125.
[11] Rogers C, Knight Jr C. An axisymmetric linear/highorder finite element for filamentwound composites—I. Formulation and algorithm. Computers & structures 1988; 29(2): 26571.
[12] Krishna Murthy A, Reddy T. A higher order theory of laminated composite cylindrical shells. 1986.
[13] Ye J, Soldatos K. Threedimensional vibration of laminated cylinders and cylindrical panels with symmetric or antisymmetric crossply layup. Composites Engineering 1994; 4(4): 42944.
[14] Kant T, Menon M. Higherorder theories for composite and sandwich cylindrical shells with C_{0} finite element. Computers & structures 1989; 33(5): 1191204.
[15] Timarci T, Soldatos K. Comparative dynamic studies for symmetric crossply circular cylindrical shells on the basis of a unified shear deformable shell theory. Journal of Sound and Vibration 1995; 187(4): 60924.
[16] J. C. Theory of thick laminated composite shallow shells. PhD Thesis, The Ohio State University 1993.
[17] Leissa A, Chang JD. Elastic deformation of thick, laminated composite shells. Composite structures 1996; 35(2): 15370.
[18] Qatu MS. Accurate equations for laminated composite deep thick shells. International Journal of Solids and Structures 1999; 36(19): 291741.
[19] Lam K, Qian W. Free vibration of symmetric angleply thick laminated composite cylindrical shells. Composites Part B: Engineering 2000; 31(4): 34554.
[20] Icardi U, Ruotolo R. Laminated shell model with secondorder expansion of the reciprocals of Lamé coefficients H_{a}, H_{b} and interlayer continuities fulfilment. Composite structures 2002; 56(3): 293313.
[21] Byon O, Nishi Y, Sato S. Optimizing lamination of hybrid thickwalled cylindrical shell under external pressure by using a genetic algorithm. Journal of Thermoplastic Composite Materials 1998; 11(5): 41728.
[22] Chandrashekhara K, Nanjunda Rao K. Approximate elasticity solution for a long and thick laminated circular cylindrical shell of revolution. International journal of solids and structures 1997; 34(11): 132741.
[23] Ganapathi M, Patel B, Pawargi D. Dynamic analysis of laminated crossply composite noncircular thick cylindrical shells using higherorder theory. International journal of solids and structures 2002; 39(24): 594562.
[24] Khare RK, Rode V, Garg AK, John SP. Higherorder closedform solutions for thick laminated sandwich shells. Journal of Sandwich Structures and Materials 2005; 7(4): 33558.
[25] Zhen W, Wanji C. A globallocal higher order theory for multilayered shells and the analysis of laminated cylindrical shell panels. Composite Structures 2008; 84(4): 35061.