Iournal of Sound and Vibration (1981) 79(3), 397-414
STRESSES IN SHELL STRUCTURES R. L. NELSON Engineering Sciences Division, Central Electricity Research Laboratories, Central Electricity Generating Board, Leatherhead KT22 7SE, England (Received 8 October 1980, and in revised form 15 May 1981)
In this paper the stresses obtained for various (thin) shell structures by using two types of doubly curved finite elements are compared with published information. One of the elements-a ring shell element-is designed to analyze axisymmetric structures such as cylinders and hyperboloids. The accuracy and convergence of this element is shown to be excellent. The other element-a quadrilateral she11 element-is designed to calculate stresses, mode shapes and frequencies of axisymmetric structures as well as sections of shell structures. The quadrilateral element is more versatile than the ring element. However, it is found that the convergence of the ring element is superior to that of the quadrilateral element. The resonant stresses of a hyperboloidal shell structure have been presented, and, as far as the author is aware, a similar investigation has not previously been reported in the literature. Whilst the primary purpose of the paper is to examine the usefulness of the two doubly curved elements in the analysis of shell structures, the examples considered in the text are described in detail to facilitate comparative structures by future investigators.
1. INTRODUCTION element computer programs are powerful tools in the analysis of vibrating structures. Programs are now available that can accurately predict the resonant frequencies and mode shapes of particularly complicated structures such as hyperboloids with “cutouts” [l] and cooling towers with column-support . Finite element programs can also be used to calculate stress distributions in such structures. However, in common with other theoretical methods, it is easier to obtain fairly accurate resonant frequencies for the structure being analyzed than it is to obtain accurate displacements, and the calculation of stresses is even more difficult. Doubly curved two dimensional finite elements are of appreciable value in the analysis of (thin) shell structures such as cooling towers, chimneys, fluid containers, etc. Quadrilateral shell elements (for which axisymmetry is not assumed) are more versatile than axisymmetric (or ring) finite elements, as the former can analyze both axisymmetric structures and selections of shell structures such as arches, dams, etc. However, ring elements are usually more accurate for a given number of degrees of freedom (d.o.f.). This paper investigates the use of ring and quadrilateral finite elements derived by Deb Nath [l] (reported in this journal previously) when used to calculate stresses in thin shell structures. The quadrilateral element is very similar to that derived by Woodman and Severn , and has been used to analyze hyperboloidal shells of revolution, conecylinder combinations and cooling towers [l, 41. The ring and quadrilateral elements are both derived from the same elasticity equations of Vlasov . The three orthogonal displacements (u, u, w) at the middle surface of the quadrilateral element are defined respectively, in terms of one-dimensional cubic interpolation polynomials  and specified nodal variables. There are 12 d.o.f. at each of the four nodes: u, au/as, au/&%,,
397 0022-460X/81/230397 + 18 $02.00/O
@ 1981 Academic
Press Inc. (London)
R. I.. NEISON
a2u/asaeo,v,..., i&/as a& (see the Appendix for notation). For the ring element a seventh order polynomial is employed to relate the displacement amplitudes at the middle surface to nodal variables. There are 12 d.o.f. at each of the two nodes: v,..., a3w/as3. U, aujas,a2ulas2,a3ulas3, The primary purpose of this paper is to investigate the accuracy of the stresses calculated by using the two shell elements, and to determine their usefulness in the analysis of shell structures. The examples considered in the text are well documented to allow their use, i.1 comparative studies, by other workers in the field. Following proof of the accuracy of both elements the resonant stresses of a hyperboloidal shell are presented, 2. STRESS-DISPLACEMENT
In this section a brief description of the derivation of the stress-displacement matrices for the doubly curved shell finite elements is given. If [(+(a, p, y)] are the stresses component at any point ((w,/3, y) in the shell, and cr, p, y are a tri-orthogonal set of right-hand curvilinear co-ordinates (see Figure l), the stress
resultants and couples distributed over a face of the shell with a constant value of (Y,for unit arc length of the middle surface, are given by
where h is the thickness of the shell wall. The stress resultants and couples [N] and [M] are depicted in Figures 1, 2 and 3. If [P] = [El, then
[PI = [Hz1LB1[Q, PH.
stress and shearing
Matrix Hz is obtained by integrating equations (1) and (2), and when terms of order h4 or smaller are neglected it is given by
PK/R i/R 0
sin C$ 1 a ---__ La3
gF i? -8r )
-p cos I$
0 v&(Rm’-p/2) 0 0
(4) at the middle surface and is given
1 a --
The matrix B relates the strains and displacements by Vlasov [S] (on page 250 of his text) as
1 L a2 -if-_-i 2 r* L as >
R. L. NELSON
and g are dependent on whether a ring or quadrilateral element is being considered. For the ring element, F = n and g = -1. For the quadrilateral element, F = (l/4) (a/[email protected]
) and g = +l. [a((~, p)] are the orthogonal displacements (or displacement amplitudes) at the middle surface of the element. In a finite element analysis it is obtained from the expression [a(%
where X is the finite element displacement 3.
PII = [xlrw, function and 8’ are the nodal variables.
This paper is concerned primarily with stresses at resonance which are likely to be of more significance than stresses of the non-resonating structure. The ring and quadrilateral shell elements have been incorporated into the programs RESAP  and SACTIL , respectively. 3.1.
The cylinder investigated was identical to that analyzed by Adelman, Catherines and Walton . The cylinder was deemed to be clamped at each end; its material properties are given in Table 1. The mode chosen by Adelman et al. , and in this present investigation, has three circumferential wave lengths (n = 3) and two meridional nodes (m = 2). This mode was investigated because quite complicated stresses are encountered near the clamped edges and the largest stresses occur in the vicinity of the edges. Some stresses exhibit steep gradients and therefore this problem constitutes a severe test of a finite element method. It should be noted that no attempt was made to use more elements at or near to areas of high stress concentrations, and that elements were ordered so that they were of equal meridional length. This enabled a more convenient investigation and appraisal of element performance. Whilst it is expected that the ring element would be more efficient than the quadrilateral element for the clamped cylinder, the study of the behaviour of the latter element is useful as it would indicate its probable efficiency when used to analyze rotationally periodic structures [lo] such as cooling towers with column-supports and hyperboloids with periodic cut-outs [l].
Properties of Adelman ‘s cylinder 0.3048 m (12 in)
0.254 x 10m3 m (0.01 in)
76.2 x 10m3 m (3 in)
Radius Mass density Young’s Poisson’s
7839.924 kg me3 (7.33 x 10m4 lb s2 in?) 206.8428 GN m-* (3 x 10’ lb in?) 0.3
IN SHELL STRUCTURES
The exact solutions for the cylinder obtained by Adelman et al.  together with solutions obtained by the program RESAP and SACTIL are presented in Figures 4 6-8. The symmetry of the mode in the circumferential direction was used to reduce size of the problem (see reference [l]) when using SACTIL. The exact solutions given as continuous lines.
the and the are
Meridional positlon, Lo (In) Figure 4. Meridional (axial) stress resultant, N, (Adelman’s cylinder).
The value of N, (Figure 4) predicted by both RESAP and SACTIL were indistinguishable from the values predicted by the exact solution. In the analyses, respectively, four ring elements and 18 quadrilateral elements were used. The quadrilateral element mesh had six elements along the meridian and three elements along the circumference. This mesh will be referred to as a (6 x 3) mesh and the equivalent mesh used for a doubly curved structure is depicted in Figure 5.
Figure 5. (6 x 3) quadrilateral shell element mesh for the analysis of a doubly curved structure with symmetry boundary conditions.
R. I.. NEIXIN
The results for the circumferential stress resultant, Np, are given in Figure 6. Point D is the maximum negative stress given by the exact solution, B that predicted by RESAP, A the value given by SACTIL and C the maximum negative value predicted by the tinite element method of Adelman et al. (This nomenclature is also used in other figures.) It I
. . . ., RESAP 4 elements;
stress resultant Na. - - -, SACTIL exact solution; - -, Adelman
(6 x 3) elements; et al. 
SAP (4 elements)
Figure 7. Resultant moment about meridional RESAP 4 elements; -, exact solution.
M,. - - -, SACTIL
(6 x 3) elements;
. . ~. .
is observed that the maximum negative stress resultants predicted by the three finite element methods (i.e., points A, B and C) are reasonably accurate when compared with the exact solution (point D). The oscillatory behaviour of the stresses predicted by all three finite element methods is due to employing too few elements in the analyses, especially near the stress boundary. When eight ring elements were used the stresses predicted by RESAP were observed to lie very close to the exact values (see Figure 6) and over most of the meridian the predicted stresses are indistinguishable from the exact solution. A large number (20) of ring elements were used, so as to be reasonably sure that convergence had occurred, and the calculated stresses were then seen to be coincident with the exact values over the whole meridian. (Note that even with this relatively large number of elements the computing time and region requirements were reasonable-about 150 s with 700 kbytes of storage in the IBM 370 Multi-Variable Task System.) A similar test was not conducted with the quadrilateral element as the computer storage required for the (6 x 3) mesh was nearing the limits imposed by the computing system, which was 750 kbytes. (It is important to report that this limitation would not have been effective if the method presented by Thomas [lo] had been used. For example six elements would have given exactly the same results as the (6 x 3) mesh. However, at the time of the investigation the routines required to implement the method were not available.) The variation in the values of M, (resultant moment about the meridional directionr along the meridian of the cylinder is given in Figure 7. The exact solution shows that a
-6 -7 -6 -Bs -10 -II D -12 : 1 -13 , \ -14 -: 1 L! 0
vP. ---. 8. Outer fibre meridional stress component, -, exact solution: - -, Adelman er ol. 191.
16 k 3) elements;
R. L. NELSON
stress boundary layer is present. The stress moment, M,, predicted by RESAP is now seen to be continuous at element boundaries. It will be noticed from Figure 7 that the moments predicted by SACTIL exhibit discontinuities. This is because the expressions for the moments contain second derivatives of the radial displacement and the quadrilateral element does not contain these derivatives as nodal variables. The maximum negative value of M, predicted by RESAP (point B) is grossly in error when compared to the value given by the exact solution (point D). (Note that the maximum stress predicted by SACTIL is almost coincident with the value for point D. However, this is considered to be a chance result and of no conseqeunce.) The oscillating behaviour of M, predicted by RESAP was observed to die out and the maximum negative stress approach D as the number of ring elements was increased. When fourteen ring elements were used, the predicted values of M, and the exact values were indistinguishable. The variation of the outer fibre circumferential stress component (cop) along the meridian is shown in Figure 8. Because of the very steep stress gradients near the ends of the cylinder, the stresses predicted by RESAP and SACTIL oscillate about the exact value. The maximum negative stresses predicted by Adelman et al. (point C), SACTIL (point A) and RESAP (point B) are in error when compared with the exact value (point D). 3.1.1. Convergence of the ring element From the results discussed above it appears that the ring element converges rapidly and accurately, and permits computing resources to be utilized very efficiently. If it can be verified that for a given number of elements the element yields accurate results, the element can be used with advantage for analyzing structures of engineering interest and to obtain reference values with which the results from other theoretical techniques can be compared. The convergence of the resonant frequencies calculated by using the ring element, for the cylinder clamped at both ends, is observed in Figure 9 for the n = 3, m = 2 mode. As this mode is associated with a relatively complicated stress distribution, it is a somewhat stringent test of the convergence behaviour of the element. It is seen, however, that the convergence is monotonic and rapid. Note that the resonant frequencies predicted with
Degrees of freedom
Figure 9. Convergence - - 0 - -, second eigenvalue.
values for the program
172 and 240 d.o.f. respectively were identical to the accuracy of the print-out (six significant places). The exact solution is 1159.1026 Hz for the first eigenvalue and there is a negligible difference of 0.002% between this value and that predicted by RESAP for the d.o.f. greater than 172. When the d.o.f. exceeds 100 the frequency has almost converged (error < 0.026% ). Thus for conditions that are not extreme the use of a large number of ring elements (say greater than 14) would yield very accurate results. To complete the study of the cylinder clamped at both ends, and to facilitate the use of these results by other workers, the variation of the orthogonal displacements U, u and w along the meridian are given in Figure 10. These displacements were obtained by using RESAP with 20 elements. (Note that the scale used for each component of displacement is different.)
-5 -6 -7
Meridional length Figure
u, u and w along meridian.
A< 8.03 in (4.797O) Figure 11. Spherical shell geometry and mesh for r oint loading Poisson’s ratio = 0.3; thickness = 0.1 in, mass density = 10 lb s’ in ~4.
5356in (3.198=‘--d_ ) Young’s
I B 2,612 in
= 10’ p.s.i.:
0.0025 - O.~XIO-~
along AB. -,
along AB (In)
(100 x 100) series solution
Yang , 0, w; 0, L’.
In the previous section it was shown that the resonant stresses predicted by the ring and quadrilateral shell elements for a cylinder compare favourably with published information. In this section the stresses obtained for a doubly curved shell (with meridional curvature) by using the quadrilateral element will be compared with well documented data published in the literature. The present analysis is for a static response as there is little or no information in the literature regarding the dynamic stresses generated in doubly curved thin shell structures. The dimensions and material properties of the spherical shell are given in Figure 11. (British units were used by previous workers [ll-131 and, therefore, for ease of comparison their system of units has been used in the figure.) The shell was assumed to be simply supported, and as the structure is doubly symmetrical, only a quadrant was analyzed. The finite element mesh employed for the analysis is also given in Figure 11. The centre of the shell was subjected to a point load of 100 lb. Smaller elements were used near the point load. The orthogonal displacements calculated by using SACTIL compare favourably with those obtained by Yang  who used a series solution method-see Figure 12. The results obtained by Dawe  were also found to be very similar and are therefore not reproduced in the figure. Table 2 compares the vertical deflection under the load as obtained by using several methods. Figures 13 and 14 show that good agreement has been obtained between the stresses predicted by SACTIL and the series solution. The finite element solutions of Yang and Dawe gave very similar results to that predicted by SACTIL and are not given in the figures. The couples calculated by SACTIL (see Figure 14) at the end of element 8 (point G) and at the beginning of element 9 (point H) are not coincident. This is because the
STRESSES IN SHELL STRUCTURES TABLE 2 Vertical deflection under load (spherical shell) Method
solution [ 121 (100X 100 term)
Gallagher [ 1 l]
couples increase rapidly in of the radial displacement The outer fibre stresses and to complete the study
magnitude near to the point load, and the second derivatives (w) are not explicitly defined as nodal variables. a, and ark have not been calculated by the previous workers these components of stress are presented in Figure 1 S.
The literature appears to be sparse with respect to information on the stresses generated at resonance in doubly curved structures, such as hyperboloids. (If resonant stresses in hyperboloids cannot be calculated accurately then stresses in resonating structures such as cooling towers subject to dynamic wind pressures cannot be accurately determined.) Hence, in this section the stress analysis of a hyperboloid is discussed in detail. The properties of the hyperboloid are given in Table 3. The hyperboloid was taken to be clamped at the end-circle of larger diameter. The n = 4, m = 2 mode was chosen for the analysis. This mode was not the automatic choice
Figure 13. Meridional and circumferential solution Yang ; 0, N, ; 0, N,.
along AB (in)
N, and NB along AB. -,
R. L. NELSON
* -I -2
Dimnce along A8 (in1 Figure 14. Resultant moments about meridional and circumferential directions M, and MB. series solution Yang [ 121; 0, M, ; 0, M,.
14 13 I2
U_ ; ---a- -; ro.
Figure 15. Outer fibres stresses cr, and ~a. --O-,
along AB (III)
Properties of hyperboloid a
Wall thickness Axial length of neck of hyperboloid from end-circle of smaller diameter Mass density Young’s modulus Poisson’s
Axial height of hyperboloid a and b are the coefficients
18.5928 m 2404.46
0.15 100.7872 of the hyperbolic
for the study as it is not associated with the lowest resonant frequency. (Note that a comprehensive table comparing the resonant frequencies calculated by using the quadrilateral and ring finite elements and other methods has been presented in reference [l]. In addition, the associated mode shapes for n = 0, . . . ,7 have been presented in reference .) However, as the lowest frequency of hyperbolic cooling towers are usually associated with this mode it was thought to be of some engineering interest to study it in preference to the other modes. The efficiency of convergence of the ring element when used for this particular problem is deduced from Figure 16, where the variation of both the resonant frequency and the
Figure 16. Convergence of resonant frequency and circumferential stress hyperboloid. --O---, Frequency; --C -, Na (multiply scale factor by - 1000).
17. Stress resultants
R. I- NELSON
circumferential stress IV0 with the effective number of d.o.f. is given. For more than 180 d.o.f. the convergence appears to be complete. The stress resultants for the hyperboloid are given in Figure 17 and the resultant moments are given in Figure 18. The stress components at the other surface of the structure are given in Figure 19 where it is seen that the meridional outer fibre stress a, is dominant. Note that. the above values were obtained by using 19 ring elements (236 d.o.f.) thus ensuring that convergence had occurred. These results can therefore be regarded as the correct values for the hyperboloid. (When the hyperboloid was analyzed by using only four ring elements the stresses calculated were almost identical to those depicted in Figures 17-19.)
?? ?? ?? ?
/ .’ /
64 x 107 64x10’ 64x10’
fibre stress components
A (6 x 3) quadrilateral element mesh (see Figure 5) was used to analyze the hyperboloid. The stress resultants, depicted in Figure 20, show a non-smooth variation due to inadequate convergence. However, if the stress distribution given in Figure 20 were to be superimposed on Figure 17 it would then be seen that reasonable accuracy has been achieved. (Note that NP in Figure 20 agrees least well with the values given in Figure 17, and that this component is small in comparison to N, and Nap.) The stress moments
for hyperboloid (SACTIL, 6 x 3 elements).
Figure 21. Stress moments for hyperboloid (SACTIL, 6 x 3 elements).
are given in Figure 21 and are seen to compare well with the values given in Figure except
for the values of M, at the ends.
Frequencies and stresses of axisymmetric doubly curved shell structures can be accurately and efficiently calculated by using the doubly curved ring finite element. For most investigations an element mesh employing between four and 10 elements can be used to give satisfactory results. A large number of ring elements (more than 14, say) can be used to give extreme accuracy when the stress gradients are not too steep, or good accuracy when the structure exhibits extremely steep gradients. (When 20 ring elements were used the computing time and region requirements were found to be still within acceptable limits.) When non-axisymmetric structures are analyzed the quadrilateral shell element can be used. The accuracy of the calculated resonant frequencies is good. Whilst acceptable accuracy was obtained when calculating the stresses in an axisymmetric doubly curved structure, very accurate values were not obtained as sufficient elements could not be employed due to restrictions imposed by the computing system. (This limitation, however, can be effectively removed by utilizing a technique which makes use of the rotational periodicity of the structure to decrease the size of the mesh [lo]. However, the computer routines necessary to apply this method to the calculation of accurate displacements, and hence stresses, had not been developed at the time of the investigation.) The stresses in a spherical shell due to a static point force were calculated and found to be in good agreement with published information. The stresses of hyperboloidal shell at resonance have been bresented and, as far as the author is aware, similar information has not previously been reported in the literature. The quadrilateral elements, validated herein, are now being applied to the calculation for resonant stresses in cooling towers (with column-supports explicitly defined). In this work the method of analysis for rotationally periodic structures described by Thomas [lo] is being utilized and results will be published in due course.
ACKNOWLEDGMENTS The work was carried out at the Central Electricity Research head) and the paper is published by permission of the Central Board.
Laboratories (LeatherElectricity Generating
REFERENCES 1. J. M. DEB NATH 1974 Journal of Sound and Vibration 33, 79-101. Free vibration, stability and “non-classical modes” of cooling tower shells. 2. T. Y. YANG, C. T. SUN, H. Lo and J. L. BOGDANOFF 1976 IEEE-ASME Joint POWU Generation Conference, Paper No. 76-JPGC-PWE- 11,2-12. Some results on dynamic characteristics of major components of 1200 MW fossil fuel steam generating plant. 3. N. J. WOODMAN and R. T. SEVERN 1970 Symposium on Structural Dynamics, Loughborough University of Technology, Paper No. A5. A double-curvature shell finite element and its use in the dynamic analysis of cooling towers and other shell structures. 4. J. M. DEB NATH 1974 Journal of Sound and Vibration 36, 253-272. Use of higher order displacement functions in the free vibration analysis of shells of revolution having meridional singularities. 5. V. Z. VLASOV 1964 NASA Publication N64- 19883 General theory of shell and its application in engineering.
STRESSES IN SHELL STRUCTURES
6. L. A. SCHMIT, F. K. BOGNER and R. L. FOX 1968 American Institute of Aeronautics and Astronautics Journal 6, 781-791. Finite deflection structural analysis using plate and shell discrete elements. 7. R. L. NELSON 1977 Central Electricity Research Laboratories, Leatherhead, Publication No. RDILfP 9177. REASP: A program to calculate resonant stresses, frequencies and mode shapes of axisymmetric structures. 8. R. L. NELSON 1979 Central Electricity Research Laboratories Leatherhead, Publication No. RD/L/P l/79. SACTIL: A program to calculate frequencies, mode shapes and stresses of shell structures. 9. H. M. ADELMAN, D. S. CATHERINES and W. C. WALTON 1970 American Institute of Aeronautics and Astronautics Journal 8, 462-468. Accuracy of nodal stress calculations by the finite element method. 10. D. L. THOMAS1979 International Journal for NumericalMethods in Engineering 14. 81-102. Dynamics of rotationally periodic structures. 11. R. H. GALLAGHER 1966 Bell Aerosystems, Buffalo Report, No. 85900-90211. The development and evaluation of matrix methods for thin shell structural analysis. 12. T. Y. YANG 1973 Journal of the Engineering Mechanics Division, Proceedings of the American Society of Civil Engineers 99, EMl, 157-181. High order rectangular shallow finite element. 13. D. J. DAWE 1974 Conference on Finite Elements for Thin Shells and Curved Members at University College, Cardiff. Some high order elements for arches and shells. 14. J. M. DEB NATH 1972 Central Electricity Research Laboratories (Leatherhead), Publication No. RD/L/N 209/72. A finite element for the analysis of the vibration and stability of shells of revolution.
matrix relating strains at reference (middle) surface of shell to displacements at reference surface Young’s modulus is n for a ring element and is (l/G) a/&? for a quadrilateral element stress-displacement matrix length of the meridian of an element couple acting in the plane through the meridian and normal to middle surface couple acting in the plane through the circumferential line and normal to the middle surface twisting moment; couples act in opposite senses in the plane along the meridian and normal to the middle surface direct stress in meridional direction direct stress in circumferential direction shear stress in the circumferential direction radius of revolution (shortest distance from the axis of revolution to the middle surface of shell; i.e., the radius of a parallel circle) tri-orthogonal set of “right-hand” curvilinear co-ordinates: cy is the direction along the meridian of the shell structure (in a downward direction), /? is the direction along the circumference of the shell structure and y is the direction through the thickness of the shell structure (zero at the middle surface of the shell) displacement at any point on the middle surface (y = 0) nodal variables angle subtended by the arc of a parallel circle =el* = R x 0 (in this paper &, is synonymous with /3) = (1 -P)/2 Poisson’s ratio stress component along meridian stress component along circumference shear stress component angle the normal to the shell wall makes with the axis of revolution angle subtended by the parallel meridional edges of the finite element at the axis of revolution
a, b f
6 m n
s J-4u, w
R. L. NELSON
coefficients of hyperbolic equation is - 1 for ring element and is + 1 for quadrilateral elemen thickness of shell wall = h3/12 meridional nodes harmonic circumferential wave number = - (r-’ -R-l sin 4) = (r-l+ R-’ sin 4) first principal radius of curvature associated with (Y (or s) distance along the meridian (it is synonymous with (Y) = s/L displacements in the (Y,@ and y directions respectively