INTRODUCTION
Thin shells and plates are very common in structural constructions. In the
design of thin shells (or plates) buckling is an important factor to be considered.
Buckling may be overall (global) or local. Buckling of axisymmetric shells,
especially cylinders, received a good amount of attention in the past (Timoshenko
and Gere, 1963; Bushnell, 1985; AlQablan,
2010). A lot of theoretical work to predict the buckling load of these structures
is available in the literature. With finite element method, it is possible to
analyze more complex shells and other structures (Wood and
Zienkiewicz, 1977; Boumechra and Kerdal, 2006; Bagchi
et al., 2007; ElKafrawy and Bagchi, 2007).
An insight on the local buckling of axisymmetric shells has been provided by
Cai et al. (2002), while Athiannan
and Palaninathan (2004) presented an experimental evaluation of stability
of shells of revolution under axial and shear stresses. Jasion
(2009) conducted linear buckling and nonlinear postbuckling analysis of
shells of revolution under external pressure and determined that convex barrel
shape shells are not stable in the postbuckling regime. Bochkarev
and Matveenko (2011) studied the dynamic and stability behavior of shells
of revolution under internal fluid pressure. Ummenhofer
and Knoedel (2000) studied the effect of boundary condition on the behavior
of steel cylindrical shells under wind load. The initial geometric imperfection
that may be present in a shell reduces the buckling capacity. There are some
experimental and analytical works on shells of revolution with geometric imperfection
available in the literature (Freskakis, 1970; Freskakis
and Morris, 1972; Frano and Forasassi, 2008). Abdullah
et al. (2008) studied the fatigue life of shells with variable amplitude
loading, while Khamlichi et al. (2010) and Bahaoui
et al. (2010) explored the effect of localized interacting defects
on the critical loads for cylindrical shells under compression. Thermal effect
of the stability of axisymmetric shells has been reported in the literature
(Bagchi and Paramasivam, 1995; Sheng
and Wang, 2010; Ghorbanpour, 2002).
In the finite element context, there are two types of buckling analysis methods normally used, linear and nonlinear. Linear buckling analysis is also called eigenvalue buckling analysis because, it gives rise to an algebraic eigenvalue problem. This method has two versions, classical and fully linearized buckling analysis. Linear buckling analysis gives good prediction of buckling load of a structure if the prebuckling rotations are negligible. However, in practice, shell structures have considerable prebuckling rotations and linear or eigenvalue buckling analysis alone is not sufficient to predict the stability limit of these structures. Full nonlinear analysis would give the exact estimate of collapse load in those cases. But full nonlinear analysis is very costly and time consuming when the size of the structural model is large.
An improved method for prediction of buckling load was suggested by Brendel
and Ramm (1980). The linear stability analysis usually performed for initial
position was repeated at different load levels on the nonlinear prebuckling
path. Thus a current estimate of the final failure load was obtained and using
this information the so called eigenvalue function was traced. They used the
characteristics of the eigenvalue function in predicting the nonlinear buckling
load. Based on this methodology Chang and Chen (1986)
proposed a scheme for predicting actual buckling load of a structure by properly
combining linear buckling analysis and geometrically nonlinear prebuckling analysis
with minimum number of load steps. They have implemented the method for beams
and general shell elements and demonstrated the effectiveness of the method
by some numerical experiments.
The simplified nonlinear buckling analysis procedure developed by Chang
and Chen (1986) for shells of revolution was extended by Bagchi
and Paramasivam (1996). The method was implemented by using conical shell
element with four degree of freedom (d.o.f) per node (Navaratna
et al., 1968) as shown in Fig. 1 which was initially
extended to solve for thermal stress and buckling capacity of shells of revolution
(Paramasivam et al., 1995; Bagchi
and Paramasivam, 1995) The basic philosophy of this method is that critical
load obtained by eigenvalue buckling analysis based on stressed structure under
a certain load level, in most of the cases, will be closer to the actual stability
limit than that obtained from an unstressed structure.

Fig. 1: 
The conical frustum shell (CFS) element 
This method can be used to determine the bifurcation point (where the load
deflection curve deviates from the primary path), if that exists. The present
article provides an overview of the stability analysis of axisymmetric shells
with implementation of the linear and simplified nonlinear buckling analysis
procedure mentioned earlier including the effect of initial geometric imperfection.
THE AXISYMMETRIC SHELL ELEMENT
A conical frustum shell element (Fig. 1) is used for the present analysis. The element has two nodal circles (i and j) and four degrees of freedom (d.o.f.) at each nodal circle. The d.o.f.s correspond to meridional displacement (u_{i} and u_{j}), circumferential displacement (v_{i} and v_{j}), normal displacement (w_{i} and w_{j}) and meridional rotation (β_{i} and β_{j}). There are three components of the field displacement, u, v and w; where, u and v are linear functions and w is a cubic function. The displacement functions are expressed in terms of the nodal d.o.f.s and Fourier harmonics, n (to represent the circumferential variation). L is the length of the element, s is the coordinate along L and ξ = s/L (a non dimensional coordinate):
The field displacements u, v and w can be expressed in terms of the element degrees of freedom using appropriate shape functions as shown in Eq. 2:
LINEAR BUCKLING ANALYSIS
In most shell structures, prebuckling stresses consist of both membrane stresses
and bending stresses and the problem basically becomes nonlinear in terms of
load and displacements. Because of this fact it is more realistic to find the
nonlinear response of the structure by step by step incremental iterative method
which is an expensive and time consuming procedure. Linear buckling analysis
is comparatively very efficient but the results obtained by this method should
be interpreted carefully.
The straindisplacement relation of a structural element in the discretized
form can be written as ε = Bδ, where B is the strain matrix. B matrix
can be broken into two, part is contributed by the linear part of the strain
and the other part is contributed by the higher order terms in the generalized
strain expression. Considering the generalized straindisplacement up to the
quadratic term B matrix can be expressed in the following form (Wood
and Zienkiewicz, 1977):
Differentiating the expression in Eq. 3 with respect to δ, we can obtain the straindisplacement relationship in the incremental form as shown in Eq. 4:
Then equilibrium equation in the incremental form can be expressed as in Eq. 5:
The solution of Eq. 4 can be approached iteratively, using methods such as, the NewtonRaphson (NR) method or its modified version. Taking appropriate variation of Eq. 6 with respect to δ we have:
where, K_{T} is the tangent stiffness matrix which can be expressed as the sum of K_{0}, the initial stiffness, K_{L}, the stiffness contribution due to geometric nonlinearlity or the second order strain as shown in Eq. 2 and K_{σ}, the stiffness due to axial stress or the geometric stiffness, as shown in Eq. 7. The expression for each of these stiffness terms are shown in Eq. 8:
In Eq. 8,
is a matrix of axial stress resultants and G represents a relation between the
radial displacements to the axial strain. D represents the constitutive matrix.
The stressstrain relationship is given by Eq. 9, where, N’s
are the inplane stress resultants, M’s are the moments and í is
the Poisson’s ratio (assuming the material to be isotropic):
The tangent stiffness matrix becomes singular (i.e., the determinant of the matrix becomes equal or very near to zero) at the critical load level, indicating that the load vector would cause the displacements to increase towards infinity. Let us assume that any arbitrary reference load is applied to the structure and it deforms to a state described by δ. In linear buckling analysis it is assumed that the incremental stiffness is a linear function of the displacements. Therefore, a positive constant λ, the structural deformation is described as λ δ for the applied load vector δP. When the critical load level is reached, the following eigenvalue equation (Eq. 10) is satisfied:
This algebraic eigenvalue problem can be solved by any standard algorithm.
The lowest eigenvalue gives the critical buckling load factor and corresponding
eigenvector describes the buckling mode. This formulation is called fully linearized
buckling analysis. If the prebuckling rotations are zero or negligible K_{L}
in Eq. 11 may be removed and the eigenvalue problem reduces
to Eq. 9 which is called classical buckling analysis:
In buckling analysis of axisymmetric shells, the stiffness matrices are formed for various a number of Fourier harmonics, n and the eigenvalues are extracted corresponding to each harmonic. The lowest eigenvalue corresponds to the critical load. Derivation of the geometric stiffness matrix for the Conical Frustum Shell (CFS) element is described in the following paragraphs.
To obtain the initial stiffness K_{0}, the small stress strain relation
(Eq. 12) is considered, while for the geometric stiffness
matrix K_{σ} and the large deflection stiffness matrix K_{L},
the strains due to large displacements (Eq. 13) are required
to be considered:
The element geometric stiffness matrices thus obtained are transformed and assembled in the global geometric stiffness matrix and the eigenvalue problem is solved for classical buckling analysis. For fully linearized buckling analysis K_{L} matrix is also included.
INITIAL GEOMETRIC IMPERFECTION
Initial imperfection in a structure can be treated as initial displacements in the structure causing some initial strains. In that case the complete strain field for a structural element can be expressed as:
where,
is the normal component of initial imperfection in the geometry of the shell
(assuming that the imperfection is predominantly in the normal direction, not
in other two directions). It is assumed that the circumferential variation is
proportional to the buckling mode shape and can be expressed in terms of Fourier
series in case the buckling mode shape is unsymmetric.
NONLINEAR BUCKLING ANALYSIS
Linear eigenvalue buckling formulation may also be applied to structures having
a nonlinear prebuckling path with a bifurcation or limit point. In this case
the lowest eigenvalue is more or less satisfactory estimate of the ultimate
load. This eigenvalue analysis can be repeated for different load levels producing
an estimate of the buckling load in an advanced position. If the function of
the lowest eigenvalue is monitored parallel to the nonlinear loaddeflection
path it gives an indication of the distance between the critical load and the
load level reached. Finally the eigenvalue function (also called the predictedload
curve (Chang and Chen, 1986) intersects the load deflection
curve in the critical point (Brendel and Ramm, 1980)
as shown in Fig. 2.
From previous discussions it is clear that the linear (classical and fully
linearized) buckling analysis alone cannot give very useful prediction about
the actual critical load of the structure. To improve the the approximation
linear buckling analysis based on a stressed structure under a certain load
level P_{base} may be performed (Chang and Chen,
1986). The buckling load so obtained in most of the cases will be closer
to the actual stability limit than that obtained from an unstressed structure.
To perform this type of analysis, the base load level, P_{base} is
first applied in increments or single step and a conventional nonlinear analysis
using NewtonRaphson method is performed. The tangent stiffness matrix, K_{T}
corresponding to that equilibrium state is formed and as a succeeding step a
trial load increment is applied and based on the incremental stresses and displacements
K_{σ} and K_{L} are formed.

Fig. 2(ab): 
Scheme for simplified nonlinear buckling analysis: (a) the
predicted load curve and (b) determination of the collapse load 
At this stage the eigenvalue problem is solved to find the proper load level
beyond the fundamental load level, P_{base} at which the tangent stiffness
matrix of the structure becomes singular. If the critical load factors are plotted
against the base load factors upon which the eigenvalue extraction are performed,
P_{base} in the predictedload curve (same as eigenvalue function on
PV curve) will always terminate at the point of coordinate (P_{cr},
P_{cr}).
The combined nonlinear and linear buckling analysis outlined above can be applied
with either the fully linearized or classical buckling analysis. First, a linear
buckling analysis on a unstressed structure is performed and an approximate
critical load, P_{cr1} is obtained; then, at certain base load level
P_{base} another linear buckling analysis is performed to obtain a new
prediction of buckling load, P_{cr2}. Finally, on a (P_{base}P_{cr})
graph (Fig. 2b), extrapolation from point (0, P_{cr1})
and (P_{base}, P_{cr2}) is made to find its intersection with
the 45° line. This approach for critical load prediction requires significantly
less computation than that for conventional geometrically nonlinear analysis.
It may provide a suitable guideline for applying load increments when a loaddeflection
curve is required to be traced. This procedure referred to as simplified or
Quasi Nonlinear Buckling Analysis (QNBA) method by Bagchi
and Paramasivam (1996).
NUMERICAL EXAMPLES
Truncated hemispherical shell subjected to axial tension: This is an
interesting problem. The geometry and loading are as shown in Fig.
3. Under the action of axial tension, compressive stress develops in the
hoop direction. This compressive hoop stress leads to buckling of the shell.
This problem was originally reported by Yao (1963). He
has presented experimental as well as theoretical results. His theoretical analysis
was based on Galarkin's approach. Yao found a considerable discrepancy between
theoretical and experimental buckling loads. The present analysis supports the
theoretical value of critical load obtained by Yao (1963).
Fully Linearized Buckling (FLB) analysis improves the result to a large extent.
If the effect of imperfection is considered the value of buckling load obtained
here is comparable to the experimental result. The values of different critical
loads and number of circumferential waves are given in Table 1.
The truncated hemispherical shell problem (Fig. 3) is solved
considering initial geometric imperfection and it is found that if the maximum
deviation of the shell surface from the perfect geometry is equal to 0.4 times
the thickness of the shell, the buckling load is 11.73x10^{4} N/m (n
= 40). This value is comparable to the experimental value of the buckling load
reported by Yao (1963) which is 10.76x10^{4}
N/m (n = 39) (Table 1). The classical buckling mode shape
is presented in Fig. 4.

Fig. 3: 
Truncated hemispherical shell under tension 

Fig. 4: 
Buckled shape of the truncated hemisphere (only a 45°
segment is shown) 

Fig. 5: 
Axially loaded cylinder 
Table 1: 
Buckling load of a truncated hemispherical shell under tension 

*Values in the parenthesis indicates value of n. The critical
load is in 10^{4} lb/in 
Cylinders under axial load: Buckling of cylinders received a lot of
attention of the engineers because of the extraordinary discrepancies between
the buckling loads obtained by theory and experiment. Buckling of cylinders
is greatly influenced by initial geometric or load imperfection and residual
stresses in it. Here, the behaviour of the perfect cylinders is discussed and
the closed form solutions are used to check the results obtained by linear buckling
analysis theory. For perfect cylinders with simply supported edges, classical
solutions are available (Timoshenko and Gere, 1963).
Finite element analysis agrees well at the above mentioned points. Results of axially compressed cylinders (Fig. 5) with different L/R and t/R ratios are presented in Fig. 6 and Table 2, respectively. From these tables it is clear that L/R ratio has very little effect (unless the cylinder is very short or very long) on critical load compared to the t/R ratio. Buckled shape of an axially loaded simply supported cylinder with modulus of elasticity E = 2x10^{5} MPa, v = 0.3 and L/R = 0.4 is shown in Fig. 7 as an example.
A spherical shell with external pressure and initial geometric imperfection:
A spherical shell with geometry as shown in Fig. 8 is considered.
E = 2x10^{5} Mpa and v = 0.3. The shell is subjected external pressure.

Fig. 6: 
Critical load for axially loaded cylinders 

Fig. 7: 
Buckled shape of a typical axially loade 

Fig. 8: 
A spherical shell with geometric imperfection 
Table 2: 
Effect of R/t ratio on the critical load of axially loaded
cylinders 

Initial geometric imperfection of the shell is considerd parallel to the buckling
mode.For different values of /t
and P_{i}/P_{0} ratios are shown in Fig. 9.
Here,
the amplitude of initial geometric imperfection with respect to the perfect
meridian of the shell, t is the thickness of the shell, P_{i} is the
buckling load for the imperfect shell and P_{0} is the buckling load
for the perfect shell. Similar types of problems were studied by Freskakis
(1970) and Freskakis and Morris (1972). Although
he assumed different shape of imperfection, those results are also presented
here as reference.
It is observed from Fig. 9 that the buckling load reduces with the magnitude of initial geometric imperfection and the rate of reduction in buckling load increases with the increase in the magnitude of imperfection.
Spherical cap subjected to axisymmetric ring load: The geometry of the
spherical cap is shown in Fig. 10. Geometrically nonlinear
analysis of this shell was presented by Wood and Zienkiewicz
(1977). The limit load is obtained from the graph presented by them (Plim
= 336.0 N) and is compared with the collapse load obtained using simplified
nonlinear buckling procedure QNBA (Bagchi and Paramasivam,
1996) which is (Pcr* = 361.3 N) and found to be significantly close (deviation
of QNBA result is about 7.5% from actual axisymmetric collapse load).

Fig. 9: 
Effect of geometric imperfection on the critical load of
the spherical shell 

Fig. 10: 
A spherical cap with ring load 

Fig. 11: 
Effect of load ring diameter on the collapse load 
It should also be noted that, axisymmetric shell subjected to axisymmetric load may not show axisymmetric collapse. The present analysis method is applied for other buckling mode with different numbers of circumferential waves and a still lower value of buckling load is recorded. For nonaxisymmetric buckling Pcr* = 272.3 N. Figure 11 shows the variation of the collapse load with the position of the ring load. As the ring diameter increases, the collapse load increases.
CONCLUSIONS
The study provides a review of buckling behaviour of shells of revolutions with a systematic development of the conical shell elements and detailed formulation including geometric imperfection. It also presents the analysis of a number of cylindrical and spherical shells and compares the experimental collapse loads with the numerical results whenever available and performs the sensitivity of length, radius or length on the buckling load as appropriate. It is found that the numerical results albeit slightly unconservative can provide a good idea about the actual collapse load and the simplified nonlinear buckling analysis produces a reliable estimate of collapse loads when prebuckling rotations are significant. However, further work is necessary to study the effectiveness of the QNBA method for axisymmetric shells with more sophisticated elements. Based on the work presented in this article, the following conclusions are made:
• 
The twonoded conical frustum shell element is found to be
very efficient for stability analysis in spite of its inherent limitations
in modelling curved meridians. In the case of a curved meridian, very small
elements should be used. Curved elements are not studied in this study 
• 
The classical buckling analysis is faster than the fully linearized buckling
analysis. But the classical buckling analysis gives good prediction of actual
buckling load only in the case of structure with membrane stresses alone
or with negligible bending stresses 
• 
The fully linearized buckling analysis gives lower buckling load than
that obtained by the classical buckling analysis. But sometimes, it may
also give very conservative result when prebuckling rotations are higher 
• 
Convergence of the critical buckling load obtained by the classical buckling
analysis with increasing number of elements is from the top and the buckling
load predicted by this method is always in the higher side and unconservative 
• 
The nonlinear buckling load prediction by the simplified nonlinear analysis
is found to be straight forward and very much economical. It gives reasonably
accurate estimate of the actual collapse load of a structure with minimum
effort 
ACKNOWLEDGEMENT
Support of the Natural Sciences and Engineering Research Council of Canada (NSERC) is greatly appreciated.