ProductsAbaqus/Standard 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 loaddeflection analysis of a model that includes imperfections. Problem descriptionThe problem consists of a long, thin, metal cylinder that is simply supported in its crosssection and loaded by a uniformly distributed compressive axial stress at its ends (Figure 1). The cylinder is sufficiently thin so that buckling occurs well below yield. (When buckling occurs in the plastic range, the problem can generally be studied numerically only by loaddeflection analysis of models that include initial imperfections. The sudden change of deformation mode at collapse causes the elasticplastic response to switch from elastic to yielding in some parts of the model and from yielding to elastic unloading at other points. Eigenvalue bifurcation predictions are then useful only as guidance for mesh design.) In the particular case studied, the cylinder length is 20.32 m (800 in), the radius is 2.54 m (100 in), and the shell thickness is 6.35 mm (0.25 in). Thus, the radius to thickness ratio for the shell is 400:1. The shell is made of an isotropic material with Young's modulus of 207 GPa (30 × 10^{6} lb/in^{2}) and Poisson's ratio of 0.3. Analysis procedureIn 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 loaddisplacement 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 predictionThe 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 crosssection with $n=$ 2 and $n=$ 3) and m longitudinal halfwaves, 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 onequarter 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 quarterwave model to represent the entire cylinder in the circumferential direction. A quarterwave circumferentially subtends an angle of $\alpha =\pi /2n.$ Likewise, we need only model onehalf 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 lowerorder 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 halfwaves m can have odd or even values, the midlength edge is alternately modeled with symmetry and antisymmetry boundary conditions. Loaddisplacement analysis on imperfect geometriesThe example is continued by performing an incremental loaddeflection 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 loaddeflection 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 loaddeflection 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. Results and discussionThe results for both analyses are discussed below. Eigenvalue buckling predictionThe analytical solution given by Timoshenko and Gere assumes that the buckling eigenmode has n lobes or waves circumferentially and m halfwaves longitudinally and provides a critical stress value for each combination of m and n. The mode that gives the minimum critical stress value will be the primary buckling mode of the shell: which mode is critical depends on the thickness, radius, and length of the cylinder. For the particular case studied here, the dependency of the critical stress values on m and n is illustrated in Figure 4: each node on the surface represents a possible buckling mode. Table 1 shows the numerical values of these critical stresses for a number of mode shapes. For this geometry the minimum critical stress corresponds to a mode shape defined by $m=$ 1 and $n=$ 4; that is, one halfwave along the cylinder and four full waves around the circumference. Figure 5 shows the (1, 4) buckling mode shape predicted with the mesh of S4R5 elements. The $m=$ 1, $n=$ 0 mode corresponds to buckling of the cylindrical shell as an Euler column: for this mode the critical stress is more than 250 times the critical stress for $m=$ 1, $n=$ 4. For small numbers of axial halfwaves (m) the critical stress changes rapidly with respect to the number of circumferential lobes (n). However, for higher values of m and n the critical stresses are not very much higher than the critical stress for $m=$ 1, $n=$ 4 and do not vary much from mode to mode, as can be observed in Figure 4 and Table 1. This behavior exhibits itself in the finite element solutions, as shown—for example—in Table 2, where the results for element type S9R5 are given and compared to the analytical results of Timoshenko and Gere. The mode numbers (values of n and m) given in that table are estimated visually from inspection of deformed configuration plots of the eigenmodes. In several cases no identification is given (the mode number is listed as ``*''), because the mesh is too coarse to define any mode. As an example, consider the mesh for $\alpha =\pi /6$, which allows for an odd number of halfwaves in the longitudinal direction. This mesh can yield eigenvectors that correspond to the $\left(m,n\right)$ mode shapes (3,1), (3,3), (3,5), $\mathrm{\dots}$ or (6,1), (6,3), (6,5), $\mathrm{\dots}$ However, as described earlier, the eigenvalues do not show an ascending pattern with the number of lobes either in the circumferential or longitudinal direction because of the geometry of this problem. Abaqus will estimate the eigenvalues in ascending order, from the closest eigenvalue to zero, unless a shift point is defined. For this case the analytical solution shows that the lowerorder modes (among those that can be represented by the mesh) have very large eigenvalues: the eigenvalues reduce steadily as the number of longitudinal halfwaves increases (see the analytical solution given in Figure 4 and Table 1), approaching a slightly higher value than the critical stress for $m=$ 1, $n=$ 4. Thus, for $\alpha =\pi /6$, the number of longitudinal halfwaves in the eigenmodes corresponding to the lowest critical stress is very large; and, since the critical stresses for all of these highorder longitudinal eigenmodes are so similar, the eigenmode is rather indeterminate. The finite element mesh, however, has a fixed number of nodes longitudinally and cannot represent these very high numbers of halfwaves with any amount of clarity. Thus, the eigenvector plots show many longitudinal modes—obviously too many for the mesh to represent accurately. It should be emphasized that these remarks apply in the context of this case only. Nevertheless, the discussion offers some useful insight into more general problems of this class and illustrates some of the difficulties that can be encountered in buckling analysis. The critical stress values in Table 2 to Table 4 for the various mode shapes correlate well with the analytical solution. Figure 6 compares the eigenvalues obtained with different shell elements with the analytical solutions. Element type S9R5 provides the most accurate results among the shell elements studied. The accuracy of this element is particularly evident in the critical stresses corresponding to the higherorder modes. S4R5 and S4R elements predict somewhat higher critical loads than S9R5. STRI3 provides stiffer solutions compared to the quadrilateral elements due to the constant membrane strain representation. The element STRI65 results correspond very closely with the analytical solutions. This element can represent linear stress variation (both in membrane and bending modes) and does not have any hourglass modes. Therefore, STRI65 is a robust and efficient element. In general, STRI65 is a good choice, particularly in problems that need very accurate modeling. A close examination of the analytical solution reveals that there are several hundred modes for which the critical stress is within 15% of the ($m=$ 1, $n=$ 4) critical stress. Therefore, this example provides a severe test of the ability of the eigenvalue algorithm to predict nearly equal eigenvalues with distinctly different eigenvectors. Loaddisplacement analysis on imperfect geometriesFigure 7 shows the applied load against the axial displacement of the node at a corner of the mesh plotted for the different initial imperfection values. The figure shows that the peak load is essentially the same as that predicted by eigenvalue analysis for the smaller initial imperfections (1% and 10% of the thickness). The larger imperfection (100% of thickness) reduces the peak load by about 12%. The analysis is completed with relative ease for an extensive portion of the postbuckling response. Figure 8 shows the deformed shape of the cylinder well into the postbuckling response. The particular case shown has an initial imperfection of 1% of the thickness. The development of the postbuckling $n=$ 4, $m=$ 1 mode is very apparent. Higher axial modes are also evident: these may be mesh dependent but are not investigated further here. Input files
References
Tables
FiguresFigure 1. Cylindrical shell with uniform axial loading.
Figure 2. Crosssection deformation corresponding to n=2
and to n=3.
Figure 3. S4R5 mesh for eigenvalue buckling prediction.
Figure 4. Critical stress for various buckling modes.
Figure 5. Buckling mode shape (m=1,
n=4).
Figure 6. Critical stress versus subtended angle of quarterwave model.
Figure 7. Applied load (normalized) versus axial displacement of an end
node.
Figure 8. Postbuckled deformation (initial imperfection of 1% of
thickness).
