Buckling of a cylindrical shell under uniform axial pressure

In general, shell buckling stability studies require two types of analysis. First, eigenvalue analysis is used to obtain estimates of the buckling loads and modes. Such studies also provide guidance in mesh design because mesh convergence studies are required to ensure that the eigenvalue estimates of the buckling load have converged: this requires that the mesh be adequate to model the buckling modes, which are usually more complex than the prebuckling deformation mode. Using a mesh and imperfections suggested by the eigenvalue analysis, the second phase of the study is the performance of load-displacement analyses, usually using the Riks method to handle possible instabilities. These analyses would typically study imperfection sensitivity by perturbing the perfect geometry with different magnitudes of imperfection in the most important buckling modes and investigating the effect on the response.

The key aspect of the eigenvalue analysis is the mesh design. For the particular problem under study we know that the critical buckling mode will be a displacement pattern with n circumferential waves (Figure 2 shows a cross-section with $n=$ 2 and $n=$ 3) and m longitudinal half-waves, and we need to determine the values of n and m that represent the lowest critical stress. One approach would be to model the whole cylinder with a very fine mesh and to assume that we can then pick up the most critical mode. This approach would be computationally expensive and is not needed in this case because of the symmetry of the initial geometry. We need to model only one-quarter of a circumferential wave: the combination of symmetry boundary conditions at one longitudinal edge of this circumferential slice and antisymmetry boundary conditions at the other longitudinal edge during the eigenvalue extraction allows this quarter-wave model to represent the entire cylinder in the circumferential direction. A quarter-wave circumferentially subtends an angle of $\alpha =\pi /2n.$ Likewise, we need only model one-half of the axial length, using either symmetry or antisymmetry at the midplane, depending on whether we are looking for even or odd modes $m.$

With this approach it is necessary to perform several analyses using different values of $\alpha$ and symmetry or antisymmetry at the midplane, instead of a single analysis with a very large model. Several small analyses are generally less expensive than one large analysis, since the computational costs rise rapidly with model size. In this particular example we consider the variation of $\alpha$ in the range of $\pi /6$ to $\pi /20$, corresponding to the range $n=$ 3 to $n=$ 10. We do not consider the cases of $n=$ 1 and $n=$ 2 because we know that these will not give lower critical loads.

The mesh chosen for the analysis of such a segment of the cylinder, using element type S4R5, is shown in Figure 3. Similar meshes with element types S4R, STRI3, STRI65, and S9R5 are also used. For the triangular elements each quadrilateral shown in Figure 3 is divided into two triangles. The meshes using element types S9R5 and STRI65 have half the number of elements in the circumferential and axial directions as the meshes using the lower-order elements. No mesh convergence studies have been done, but all the meshes and elements used give reasonably accurate predictions of the critical load.

Eigenvalue buckling analysis is performed with Abaqus by first storing the stiffness matrix at the state corresponding to the “base state” loading on the structure, then applying a small perturbation of “live” load. The initial stress matrix resulting from the live load is calculated, and then an eigenvalue calculation is performed to determine a multiplier to the “live” load at which the structure reaches instability. In this example there is no load prior to the “live” load; therefore, the eigenvalue buckling (Eigenvalue buckling prediction) is the only step. During the buckling procedure one longitudinal edge has symmetry boundary conditions, and the other has antisymmetry boundary conditions, as shown in Figure 3. With these constraints a mesh subtending an angle of $\alpha$ can model modes with $n=\pi /2\alpha ,\pi /\alpha ,2\pi /\alpha \mathrm{\dots }$ waves around the circumference of the cylinder. However, during the calculation of the initial stress matrix, both longitudinal edges must have symmetric boundary conditions (because the prebuckling response that creates this stress stiffness is symmetric). Thus, the boundary conditions associated with the “live” loading are specified under one boundary condition, and the boundary conditions associated with the buckling deformation are defined under a second boundary condition. If the second definition is not given, the boundary conditions are the same for the loading and for the buckling mode calculation. The loaded edge is simply supported. Since the number of longitudinal half-waves m can have odd or even values, the midlength edge is alternately modeled with symmetry and antisymmetry boundary conditions.

The example is continued by performing an incremental load-deflection analysis using the modified Riks method. For some problems linear eigenvalue analysis is sufficient for design evaluation, but if there is concern about material nonlinearity, geometric nonlinearity prior to buckling, or unstable postbuckling response (with associated imperfection sensitivity), the analyst generally must perform a load-deflection analysis to investigate the problem further.

The mesh used for this phase of the analysis consists of eight rows of elements of type S4R5 in the circumferential direction between symmetry lines. (In the eigenvalue analysis antisymmetry boundary conditions are used, since the analysis is a linear perturbation method. But this load-deflection study allows fully nonlinear response, so the antisymmetry assumption is no longer correct.) Twenty elements are used along the length of the cylinder.

An imperfection in the form of the critical buckling mode (obtained in the previous analyses of the example) is assumed to be the most critical. The mesh is, therefore, perturbed in the radial direction by that eigenmode, scaled so that the largest perturbation is a fraction of the shell thickness. The studies reported here use perturbations of 1%, 10%, and 100% of the thickness. The following examples demonstrate two methods of introducing the imperfection.

The first method makes use of the model antisymmetry and defines the imperfection by means of a Fortran routine that is used to generate the perturbed mesh, using the data stored on the results file written during the eigenvalue buckling analysis. bucklecylshell_stri3_n4.inp shows the input data for the buckling prediction, bucklecylshell_progpert.f shows the Fortran routine used to generate the nodal coordinates of the perturbed mesh, and bucklecylshell_postbucklpert.inp shows the input data for the postbuckling analysis. The meshes for the buckling prediction analysis and the postbuckling analysis are different and are described in the “Input Files” section. The postbuckling analysis is performed using the static RIKS procedure (Unstable collapse and postbuckling analysis).

The second method uses a geometric imperfection to define the imperfection, which requires that the model definitions for the buckling prediction analysis and the postbuckling analysis be identical. bucklecylshell_s4r5_n1.inp shows the input data for the buckling prediction, and bucklecylshell_postbucklimperf.inp shows the input data for the postbuckling analysis.