1.2.3 Buckling of a cylindrical shell under uniform axial pressure

Product: ABAQUS/Standard  

This example illustrates the use of ABAQUS to predict the elastic buckling instability of a “stiff” structure (a structure that exhibits only small, elastic deformations prior to buckling). The example is a classic case of this type of problem; a detailed analytical discussion of the problem is available in Timoshenko and Gere (1961). This analytical solution allows the example to be used for verification of the numerical results.

The structural analyst often encounters problems involving stability assessment, especially in the design of efficient shell structures. Since the shell is usually designed to carry the loading primarily as a membrane, its initial response is stiff; that is, it undergoes very little deformation. If the membrane state created by the external loading is compressive, there is a possibility that the membrane equilibrium state will become unstable and the structure will buckle. Since the shell is thin, its bending response is much less stiff than its membrane response. Such buckling will result in very large deflections of the shell (even though the postbuckling response may be mathematically stable in the sense that the structure's stiffness remains positive). In many cases the postbuckled stiffness is not positive; in such cases the collapse load generally will depend strongly on imperfections in the original geometry; that is, the structure is “imperfection sensitive.” In some cases the buckling may be only a local effect in the overall response: the shell may subsequently become stiffer again and reach higher load levels usefully with respect to its design objective. Sometimes there are many collapse modes into which the shell may buckle. For all of these reasons shell collapse analysis is challenging. This example illustrates the standard numerical approach to such problems: eigenvalue estimation of bifurcation loads and modes, followed by load-deflection analysis of a model that includes imperfections.

Problem description

Analysis procedure

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.

Eigenvalue buckling prediction

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 1.2.3–2 shows a cross-section with  2 and  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 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

With this approach it is necessary to perform several analyses using different values of 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 in the range of to , corresponding to the range  3 to  10. We do not consider the cases of  1 and  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 1.2.3–3. Similar meshes with element types S4R, STRI3, STRI65, and S9R5 are also used. For the triangular elements each quadrilateral shown in Figure 1.2.3–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, *BUCKLE (Eigenvalue buckling prediction, Section 6.2.3 of the ABAQUS Analysis User's Manual) 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 1.2.3–3. With these constraints a mesh subtending an angle of can model modes with 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 the *BOUNDARY definition, and the boundary conditions associated with the buckling deformation are defined under *BOUNDARY, LOAD CASE=2. 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.

Load-displacement analysis on imperfect geometries

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, Section 6.2.4 of the ABAQUS Analysis User's Manual).

The second method uses the *IMPERFECTION option to define the imperfection. This option 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.

Results and discussion

Input files

Reference

Tables

Table 1.2.3–1 Critical stresses versus mode shape, stresses given in GPa (from Timoshenko and Gere, 1961).

12345
075.0864.2951.8640.8132.04
1116.724.4527.8426.2523.05
21.4784.7417.8329.76910.53
30.3881.2512.3893.4784.331
40.2810.4780.9131.4171.908
50.4790.2980.4490.6810.942
694.650.3290.3090.4010.533
71.7570.4950.3140.3080.360
83.0220.7940.4140.3160.305
94.8751.2510.5100.3940.322
107.4731.8980.8780.5370.395
 
678910
025.3720.3616.5913.7111.48
119.6816.6514.1011.9910.27
210.479.9419.1908.3767.577
34.8865.1655.2285.1364.945
42.3282.6542.8783.0103.064
51.1971.4301.6251.7781.888
60.6800.8270.9661.0891.191
70.4370.5250.6160.7020.782
80.3320.3770.4310.4870.544
90.3050.3150.3390.3720.407
100.3330.3100.3050.3180.336

Table 1.2.3–2 Critical stresses – S9R5 element, stresses given in GPa.

Boundary condition at midlength of cylinder
SYMMASYMM
0.316(*, *)0.318(*, *)
0.281(4, 1)0.317(4, *)
0.316(*, *)0.299(*, 2)
0.310(6, 3)0.316(6, *)
0.315(7, 3)0.309(7, 4)
0.306(8, 5)0.316(8, 4)
0.316(9, 7)0.306(9, 6)
0.310(10, 7)0.309(10, 8)

Table 1.2.3–3 Critical stresses – S4R5, S4R elements, stresses given in GPa.

Boundary condition at midlength of cylinder
SYMMASYMM
0.327(*, *)0.327(*, *)
0.290(4, 1)0.326(4, *)
0.326(*, *)0.308(*, 2)
0.320(6, 3)0.326(6, *)
0.327(7, 3)0.319(7, 4)
0.317(8, 5)0.326(8, 4)
0.326(9, 7)0.317(9, 6)
0.322(10, 7)0.320(10, 8)

Table 1.2.3–4 Critical stresses – STRI3 element, stresses given in GPa.

Boundary condition at midlength of cylinder
SYMMASYMM
0.359(*, *)0.355(*, *)
0.285(4, 1)0.357(4, *)
0.359(*, *)0.308(*, 2)
0.321(6, 3)0.334(6, 2)
0.322(7, 3)0.321(7, 4)
0.319(8, 5)0.325(8, 4)
0.332(9, 5)0.319(9, 6)
0.324(10, 7)0.326(10, 8)

Table 1.2.3–5 Critical stresses – STRI65 element, stresses given in GPa.

Boundary condition at midlength of cylinder
SYMMASYMM
0.319(*, *)0.308(*, *)
0.280(4, 1)0.315(4, *)
0.326(*, *)0.298(*, 2)
0.309(6, 3)0.328(6, 2)
0.314(7, 3)0.308(7, 4)
0.305(8, 5)0.315(8, 4)
0.315(9, 5)0.305(9, 6)
0.309(10, 7)0.308(10, 8)


Figures

Figure 1.2.3–1 Cylindrical shell with uniform axial loading.

Figure 1.2.3–2 Cross-section deformation corresponding to n=2 and to n=3.

Figure 1.2.3–3 S4R5 mesh for eigenvalue buckling prediction.

Figure 1.2.3–4 Critical stress for various buckling modes.

Figure 1.2.3–5 Buckling mode shape (m=1, n=4).

Figure 1.2.3–6 Critical stress versus subtended angle of quarter-wave model.

Figure 1.2.3–7 Applied load (normalized) versus axial displacement of an end node.

Figure 1.2.3–8 Postbuckled deformation (initial imperfection of 1% of thickness).