Post-buckling analyses of variable-stiffness composite cylinders in axial compression

Variable-stiffness shells are thin composite structures in which the reinforcement direction is a function of its surface co-ordinates. This paper presents a numerical investigation into the buckling and post-buckling of two variable-stiffness cylinders under axial compression. Both shell walls are made from unidirectional carbon fiber slit tapes that are steered to give them a piecewise-continuous fiber-angle variation around the circumference. Dynamic analyses of the complete loading and unloading cycles are computed using a time-integrated finite element model (Abaqus). The numerical results generated herein are compared with test data and are found to be in good agreement, in terms of axial force versus end-shortening and global displacement fields. The analyses provide significant new insight into the mechanisms underpinning collapse behavior of the shells. For example, the development of the initial nonlinear buckle, its dynamic snap-through, and the formation of a post-buckled configuration are clearly visible. One effect elucidated by this investigation is the symmetry-breaking mechanism of the circumferential stiffness variation. In contrast to a constant-stiffness cylinder, in which the total strain energy is invariant to the translation of a dimple of fixed dimensions, the present structures exhibit a single dominant post-buckling mode that are associated with the formation of ‘trapped’ surface dimples. In one case, this dominant mode is found to be stable over a significant amount of further end shortening.


Introduction
Cylindrical shells are found in numerous aerospace, marine and energy structures, such as solid rocket boosters (SRBs), launch vehicle fuel tanks, and supports in offshore structures. Their ability to carry high levels of axial compression is used in many roles, where most of the structure is loaded in a pure membrane state and its efficiency is derived from the absence of through-thickness stress gradients.
Cylinders with a sufficiently thin shell wall fail due to local buckling, and in this regard they have two well-known characteristics. The first is that the development of significant radial displacements is accompanied by the rapid development of global buckling phenomena and a dynamic collapse.
Examples of this effect in the case of a constant-stiffness composite cylinders are found in Refs. [1][2][3].
The second is that, for reasons connected with the first (see ref. [4] for a detailed discussion), prediction of the loads levels at which collapse occurs is difficult when compared to well-behaved structures (e.g., Euler-Bernoulli beams or Kirchhoff-Love plates). Early attempts at deriving shell critical buckling loads were based on a pure membrane prebuckling stress state and linearized equations of neutral equilibrium. The results were elegant and compact, but failed to represent the observed behavior adequately (see the reviews by Gerard & Becker [5] and Simitses [6]). Using both experimentation and theoretical reasoning, the difficulties in predicting the buckling loads were identified and are due to (1) the presence of a prebuckling boundary layer [7], (2) a sensitivity to geometric imperfections [8,9], and (3) the presence of imperfections in the loading [10]. An important theoretical step in the nonlinear analysis of isotropic shells was made by von Kármán & Tsien [8].
Using a diamond radial displacement pattern and a compatible stress function, static equilibrium configurations were calculated in the post-buckling regime. These solutions uncovered a highly unstable response, with increasing radial displacements resulting in decreasing axial force and end shortening. In addition, a post-buckled configuration having equal end shortening to the critical value, but with smaller axial force, was found. This latter finding agreed with data from displacementcontrolled experiments in which a significant drop in load was exhibited at the point of buckling.
As the precision of experimental techniques and the power of numerical tools increased, it was found that accurate predictions of collapse behavior may be obtained as long as sufficient data of the as-built structure could be included in the model [11,12]. The main question that remained, however, was that, even though the physical behavior of a fabricated structure may be modeled accurately, how should the design proceed if the critical failure mode is so dependent on manufacturing subtleties? In both experiments and simulations a single dimple in the shell-wall was found to form where local increases in axial stress were present [13,14]. Such observations led to the conjectures of Zhu, et al. [15] and Calladine [16] that the relationship between the critical buckling stress of an isotropic cylinder and its wall thickness is closely related to the theoretical stress required to invert a dimple in a doubly curved shell. While not mathematically rigorous, their explanations fitted the trends in experimental data well and were supported by the findings of Deml & Wunderlich [17] who found that, in response to the 'worst' type of imperfection, an isotropic shell's deflected shape is a pronounced dimple. Later, Hühne, et al. [1] proposed that a robustly designed shell should be one which has the capacity to resist the applied axial load, while simultaneously resisting the action of a singular perturbation load (i.e., a dimpling force acting normal to the shell wall). A good estimate of the collapse load was found to be the one in which its sensitivity to the perturbation load is reduced by a significant amount.
Recently, buckling tests have been conducted on a relatively new class of composite structure [18].
Variable-stiffness shells are thin-walled structures in which the fiber reinforcements follow curved paths over their surfaces. The shells examined in ref. [18] had a piecewise-continuous variation in fiber angle around the circumference. The experiments yielded three points of interest. The first was that a relatively good estimate of the buckling loads (by shell expectations) was provided by a linear eigenvalue analysis (7.7% and 12.4% lower). The second was that in both shells isolated dimples were present in the post-buckled configurations immediately after buckling. The third was that, in one of the structures, a significant region of secondary stability occurred which is not a typical feature in cylinder post-buckling. This paper presents the results of a nonlinear finite element analysis counterpart study investigating novel aspects of the experimental observations. Dynamic numerical analyses are used to provide insight into the present behavior, explaining effects which are not readily observed due to the high buckling velocities of the shell wall occurring in experiments.
The specimens and experimental results are presented first, followed by a detailed description of the models developed for these analyses. Results generated using both linear and geometrically nonlinear analyses are presented for the structures with and without measured geometric imperfections. In the dynamic analyses, the complete loading and unloading paths are calculated and are in agreement with the experimental results in terms of both the deformed shape, and the measured axial loads and end displacements. The evolution of the initial collapse is presented and the origin of the instability, in both cases, is identified as the rapid growth of a single inward half-wave of the initial buckling mode. Finally, the implications of the findings on future designs are discussed.

Specimens
Two cylindrical variable stiffness shells were manufactured via the automated fiber placement (AFP) process. Details of their manufacture and design may be found in references [19] and [20], respectively. In this section a short summary of their geometry and materials is presented, and the details of the test set-up are described.
Positions on shell surfaces are defined in terms of an axial coordinate, x, and a circumferential arclength, y. We also use a Cartesian system (X,Y,Z) to define positions in three-dimensional space. That is, the variable fiber-angle is 10° at the shells' crown and keel, and transitions to 45° at the shells' sides (see Fig. 1a for the circumferential locations). In the regions between these points the fibers have constant curvature in x-y space. The resulting variable stacking sequence produces laminates with a relatively high axial bending and membrane stiffness at the crown and keel (due to the fiber's near-axial orientation) and high twist and in-plane shear stiffness at the sides. A degree of bend-twist coupling (i.e. D 16 and D 26 ) is present at all points on both shell surfaces. This effect arises due to a combination of the relatively small number of plies and the high degree of layer-wise orthotropy.
Courses computed by the AFP CAD software are shown in Fig. 2, plotted on the nominal surface geometry. A well-known characteristic of the AFP process is that as material is steered away from a straight line its effective width (in the direction perpendicular to the starting direction) varies and, as a   [21]. This effect leads to either resin-rich pockets or a thickness build-up, in turn leading to local discontinuities in the stiffness and strength of the laminate.
Two specimens, one with overlaps (Shell A) and one without (Shell B) were produced [18]. In the former case, a considerable build-up of the section resulted at the crown and keel regions, in turn resulting in a nominal as-built stacking sequence: [±45,10 3 ,-10 3 ] s there. In the latter shell, tow overlaps were avoided by dropping fiber-tows at calculated regions of overlap. Both specimens were constructed using HexPly® IM7/8552 prepreg, which has a nominal 35% volume fraction and a cured ply thickness t = 0.131 mm [22].
Geometric imperfections and laminate thicknesses were determined experimentally with a coordinate measurement machine (CMM) [19]. In the case of Shell B, the laminate ply count is approximately constant at all points on its surface (except where gaps have occurred due to dropping tows). Dividing the laminate thickness data by the constant ply count results in a mean ply thickness of 0.124 mm with a standard deviation of 0.008 mm. In the case of Shell A, the ply count is not constant over the shells surface. In order to compute the average ply thickness, the laminate thickness data was combined with the known locations of ply overlap. The mean and standard deviations of the ply thicknesses results for regions of 8, 10 and 12 plies were averaged. The average ply thickness was determined to be 0.128 mm with an average standard deviation equal to 0.006 mm. A summary of the measured shell dimensions is given in Table 1.
Structural testing of the specimens is documented in ref. [18]. Figure 3 shows the experimental setup. Axial compression was applied to the specimens with a 1300 kN test machine. The specimens were mounted between two flat platens and subjected to a controlled end shortening using a bilinear ramp profile. The maximum end shortening for Shells A and B was 4.064 and 3.048 mm, respectively. The multiple displacement transducers were required to measure relative rotations of the platens about the Y-and Z-axes (see Fig. 1a for the coordinate definitions).

Finite-element models
The finite-element package Abaqus 6.12 [23] was used to generate and solve models of the buckling and post-buckling behavior of the two shells. Meshes of `S4R' elements were produced based on both a perfect cylindrical geometry and the specimens' actual surface topologies. Data for the spatial distributions of total laminate thickness, stacking sequence, shell offsets, and the measured reference-surface geometry were obtained from the design [20], manufacturing [19] and experimental [18] studies for the shells.

Element stiffness
In order to represent a variable-stiffness shell structure with commercial shell elements, it was required that, as an approximation, each element had its own individual stacking sequence. This gave the distribution of stiffness a piecewise-constant 'patchwork' form as opposed to its real smooth and continuous form. Consequently, a relatively large number of elements was required, not only to compute an accurate approximation of the displacement fields, but to provide a good approximation of the section stiffnesses. The criterion for selecting the mesh density was convergence of the first buckling eigenvalue. This resulted in an element size of 6.35 × 6.37 mm 2 and a total number of degrees of freedom of approximately 128,000. Figure 4 shows the mesh used in the analyses.
As a consequence of the manufacturing process, the effective ply thickness (defined as the measured laminate thickness divided by the local ply count) is variable over the surfaces of both shells. Histograms of the laminate thickness measurements for each shell are given in ref. [19]. In these analyses, a uniform ply thickness, denoted t*, was used with classical laminate theory to calculate the laminate stiffness constants. The selection of this thickness value is discussed in section 4.1.
Elastic constants for the cured prepreg material were determined via the ASTM standard D695 compression test [18]. The tests were performed on flat laminates produced by the same AFP process that was used to produce the test specimens. The elastic constants E 11 , E 22 , and G 12 used in the analyses are 8% larger than the values reported in ref. [18], an increase which accounts for volume fraction differences in the coupon and shell plies (see Table 2).

Damping
The Hilber-Hughes-Taylor (HHT-α) [24] time-integration algorithm, implemented in Abaqus 6.12, was used to simulate the shell collapse behavior. The method allows artificial damping to be applied to the system while maintaining second-order accuracy over the time step.
In the experiment, structural damping is, potentially, derived from several sources, including materials (mainly the composite matrix), air-structure interaction, the test machine and boundary conditions (epoxy potting material). The time-step size in a dynamic analysis is inversely proportional to the natural vibrational frequency of the structure. The present specimens are constructed out of high specific-modulus materials and, as such, have relatively high natural frequencies. Numerical damping is applied to the HHT-α analyses in order to attenuate any oscillations, circumventing the need, difficult as it would be, to determine and represent the real sources of damping in the model.
After a number of trials, it was determined that the highest possible damping (α = -1/3) was required to reduce the number of time increments to a reasonable amount, or even in some cases achieve convergence of every time increment.

Boundary conditions
In the tests, the ends of the shell were held between two stiff platens, one of which translated while the other remained static. In the models, this attachment was represented by connecting the boundary nodes at each end (i.e., those located at X = 0 and 889 mm) to an axial control node with rigid-beam multi-point constraints (MPCs). This produced a 'wheel' effect at each end where the rigid beams acted as the 'spokes' and the control point acted as the 'hub'. At the end X = 0 the control node was fully restrained, while at the end X = L the control node was restrained from translating in the Y-Z plane. Rotational constraints on the latter control node varied depending on the analysis type.
In addition to the direct contact with the test machine, elastic constraint was imposed on the shell wall by the epoxy encasement (the affected mesh regions are shown in red in Fig. 4). The most highly Starnes [11] who used a plane strain finite-element model to calculate the local increase in a laminate's axial elastic modulus due to encasement. They found that, for a given encasement material, the increase in local stiffness depends on the elastic modulus of the laminate.
In order to determine the local axial stiffness increase associated with the shell's encasement, and hence the increase in the global axial stiffness, a similar model to Hilburger & Starnes' was produced (see Fig. 5a). The model consisted of an assembly of laminate, encasement and aluminum hoop; all in plane strain. The laminate region had length 50.8 mm in the x-direction and was assumed to be perfectly bonded to the encasement. This assembly was, in turn, assumed to be perfectly bonded to the aluminum hoop region. The mid-plane of the laminate was assumed to be a line of symmetry and the platen face was assumed to be perfectly rigid and frictionless.
The local increase in the axial stiffness of the laminate due to the encasement (per unit length) was calculated in the following way. A unit force, F, in the negative x-direction is applied at the top of the laminate. This results in compression of both the free top laminate section and the lower encased section. Assuming that the un-encased and the encased parts of the laminate behave as springs in series, their axial stiffness ratio (i.e., k encased. /k free ) was calculated as: where d 1 and d 2 are nodal displacements (shown in Fig. 5b), and η is designated the encasement stiffness ratio. The potting resin was assumed to have elastic constants given in Table 2 and the following through-thickness laminate properties were assumed to be independent of stacking sequence: E z = 9 GPa and G xz = 4 GPa. The ply thickness was assumed to be 0.128 mm.
Given that some of the fiber angles were spatially varying, it was necessary to perform the analysis for a range of stacking sequences. Three extreme stacking sequences were chosen, based on their effective axial modulus. They were those that corresponded with: the side regions of both shells; the crown/keel of Shell A; and the crown/keel of Shell B.
Results, including the calculated laminate elastic constants (calculated using classical laminate analysis) and the encasement stiffness ratios, η, are shown in Table 3. Once the local stiffness ratios had been calculated for the three extreme stacking sequences, the increase in the stiffness of both shells over their axial length (at fixed circumferential locations) was calculated. This was done by considering the total encased shell wall (50.8 mm long) and the free shell wall (838.2 mm long) to be springs in series.

Results and discussion
Numerically predicted results for the two specimens are presented in this section.

Encasement and ply thickness sensitivity
In this section the results of sensitivity analysis of the dynamic buckling predictions on variations in the encasement stiffness ratio η and the uniform ply thickness t* are presented. This investigation was conducted in order to quantify the effects of these parameters on the non-linear analysis results.
First, end-shortening was applied to each of the shells equal to 1.5 times the critical value, over a period of 600 s. The analyses were performed with η equal to small (1), intermediate (2) and large (4) values (achieved by reducing the encased elements length by a factor 1/η). This relatively large range of η was chosen to produce a correspondingly large variation in the dependent variables (i.e. the prebuckling stiffness, the buckling load and the load-drop). To simplify the analysis, the increase in stiffness was applied uniformly around the circumference. It is noted that while this boundary condition does not represent the test exactly, it was thought to be sufficiently close for the purposes of a sensitivity analysis and that its variation could only increase the influence of the encasement.  Figure 6 shows the load vs. end-shortening results for both shells with t* set equal to the mean measured values (see Table 1) and with imperfect cylindrical geometries. In both shells, the largest sensitivity to variations in η was exhibited by the prebuckling stiffness (as much as 3.9% variation for shells A and B with η in the range 1 ≤ η ≤ 4). The lowest sensitivity was exhibited by the buckling loads, the variations of which were less than one percent for both shells. This is in agreement with the classical result that the local buckling of an axially compressed cylinder is invariant of length for cylinders above a threshold value and below a length where overall Euler buckling is induced. In the case of Shell A, the drop in load immediately after buckling ranged from 59.13% (with η = 1) to 57.00% (with η = 2), a difference of 2.13%.  and a noticeable difference between the mean and the modal value. In the manufacture of Shell B tows were dropped at locations of overlap, which may account for the skewing of t he thickness distribution. Figure 7b shows the load vs. end-shortening response for Shell B with the uniform thickness set equal to t* = 0.128 mm (which is the mean value calculated for the shell without tow drops). Increasing the ply thickness by 3.23% resulted in an increase of 3.34% in the pre-buckling stiffness and a 7.15% increase in the buckling load. Similar results were obtained for Shell A, with a 3.13% increase in the uniform ply thickness (i.e. setting t* = 0.132 mm) resulting in a 3.10% increase in the pre-buckling stiffness and a 6.15% increase in the buckling load.
In the proceeding analyses, the homogeneous ply thickness a ply thickness of t* = 0.128 mm and the intermediate stiffness ratio η = 2 has been used for both shells.

Linear solutions
In ref. [18] linear buckling eigenvalues were determined for the shells based on perfect and imperfect cylindrical geometries. In addition to axial translation, data from the platen displacement transducers informed that a small degree of rotation about the Y-and Z-axes (designated ϕ Y and ϕ Z respectively) occurred in the prebuckling range. This behavior indicated that the axial load was slightly misaligned with the shell's axes of rotation. In this section, the effect of including applied bending moments in the linear buckling predictions is investigated. Slopes calculated from the experimental data were used to incorporate nodal rotations into the load case of a buckling eigenvalue analysis. The rotations were applied to the models at the control nodes located at X = L. These errors are smaller than those reported in ref. [18] because in the present work the ply material  properties have been adjusted to account for discrepancies between coupon and shell ply thicknesses.
In general, inclusion of the platen rotations reduced the buckling loads and concentrated the radial displacements at the crown. Interestingly, in the case of Shell A, the largest reduction in the predicted buckling load (2.3% compared with the pure end-shortening case) was found when including rotations about the Z-axis alone. Combining both rotations for this shell led to a 1.75% reduction in the predicted buckling load. In the case of Shell B, combined end-shortening and rotation about the Z-axis resulted in a reduction of the buckling load prediction of 2.51% (relative to a pure end shortening).

Loading and unloading cycles
Dynamic buckling simulations are now discussed. A controlled end shortening was applied to the models which ramped from zero to its maximum value (4.064 and 3.048 mm for Shells A and B respectively) over a period of 600 s. This motion was then reversed in a second step over the same time period.
Platen rotations were not included in the load case for the dynamic simulations due to the loss of proportionality of the forces and rotations in the measurements which occur in the non-linear range (see Fig. 8). Platen rotations about the Z-and Y-axes were constrained over the entire loading cycle. A separate set of analyses in which the platen rotations were unconstrained was also performed. The results of these analyses, however, contained spurious solutions and did not correlate well with the experimental results.
In the proceeding sections, direct comparisons between experimental and simulated results are made assuming models with the imperfect cylindrical geometry. Figure 10 shows a comparison of the force versus end-shortening results obtained from the experiment and the simulation. Figure 11 shows qualitative color-maps of radial displacements (experimental) and displacement magnitudes (simulation), plotted on the deformed shell geometries at different points on the loading cycle. Experimental displacement data which was used to generate these plots were obtained from digital image correlation (DIC) measurements [18].

Shell A
The buckling load predicted by the non-linear analysis was 171.75 kN which is 0.52% smaller than the measured buckling load and 0.14% smaller than the buckling eigenvalue. Figure 11a  inward dimple at the crown/keel (centered at approximately X = L/2) (see Fig. 11b). The predicted and measured load drops were 100.30 and 95.55 kN, respectively (see point B, Fig. 10b). This event is referred to herein as the first localization. Further end shortening, beyond point B, resulted in an increase in the amplitude of the first localization, a small drop in load and the development of a smaller dimple near the boundary (point C). At the end of the loading step (point D, Fig. 10b) the predicted and measured axial forces were in good agreement, as well were the surface displacement patterns.
During the unloading step, the axial force gradually increased as the localized dimples inverted and released strain energy. Towards the end of the unloading step the dimple which formed at the first localization inverted and the prebuckling part of the path was rejoined. In the experimental force Figure 11. Radial displacement (experimental) and displacement magnitude (simulation) color-maps of Shell A at different points on the loading cycle. The color-code for the experimental data is: red = +ve (outward), green = 0 and violet = -ve (inward).
versus end-shortening results a gradient equal to that of pre-buckling is adopted but with a shift in end-shortening. This hysteretic effect was possibly due to backlash in the test machine and/or permanent deformation of the specimen. Figure 12 shows a comparison of the force versus end-shortening results obtained from the experiment and the simulation. Figure 13 shows qualitative color maps of radial displacements (experiment) and displacement magnitudes (simulation), plotted on the deformed shell geometries at different points on the loading cycle.

Shell B
The buckling load (point A, Fig. 12b) predicted by the nonlinear analysis is 71.22 kN, which is 6.85% smaller than the measured buckling load and 1.54% larger than the buckling eigenvalue. Figure   13a shows a comparison of the experimental and simulated radial displacements. The radial displacements in both cases increase in amplitude near the shell ends. In contrast to the simulated results, however, the experimental displacements were larger at the end X = 0 than at X = L. This asymmetry may be due to platen rotations which were observed during the test (see Fig. 8). The numbers of axial half-waves in the experiment and in the simulation were 26 and 23, respectively.
The peak load was followed by a drop in load of 20.60 kN (experimental) and 20.03 kN (simulated). In contrast to the behavior of Shell A, the occurrence of first localization (point B, Fig.   12b) was accompanied by a region of significant secondary stability (this is discussed further in At approximately 150% of the critical end shortening a secondary collapse occurred. This was accompanied by the formation of two further localized dimples. Good agreement is noted in Fig. 12 between the measured and predicted load levels. After this event, the good overall agreement between experiment and simulation ended. At point C the predicted radial displacement pattern is qualitatively similar to that observed in the experiment. The prediction of the axial load at this point, however, is quite erroneous (-27.38% error). The model based on the perfect cylindrical geometry failed to predict this tertiary state (point C, Fig. 12b) and instead adopted a configuration with 12 circumferential Figure 13. Radial displacement (experimental) and displacement magnitude (simulation) color-maps of Shell B at different points on the loading cycle. The color-code for the experimental data is: red = +ve (outward), green = 0 and violet = -ve.(inward).
waves (point E, Fig. 12a). This departure in modeling accuracy is thought to reflect the model's idealization whereby a constant thickness is used instead of the actual periodic thickness variation. At increased levels of end shortening, the model with an imperfect geometry exhibited a third collapse (which was not evident in the experiment) also with a large number of circumferential waves (point D, Fig. 12b).
During the unloading step, the load gradually increases until, at a value of end shortening approximately equal to the critical one, the same configuration as the first localization is adopted (lower point B, Fig. 12b). The structure then unloads along a path with equal gradient to that of the prebuckling stage.

Initial collapse process
One of the advantages of performing a dynamic simulation is that displacement fields which would be changing rapidly may be examined in detail by viewing the results frame-by-frame. In this section results of the development of the initial collapse are presented. Figure 14 shows the incremental changes in Shell A's deformed geometry occurring immediately after the buckling load has been reached (i.e. at points on the line A-B of Fig. 10b). At the buckling load (point A, Fig. 10b) a harmonic radial displacement field is observed with a maximum radial displacement of approximately 0.68 mm (w/h = 0.34). This deformation is of the same order of magnitude as the displacements measured in the test (0.25 mm, w/h = 0.12). Interestingly, the buckling behavior exhibited by this structure is akin to that of a rod on elastic foundation [25]; the main characteristics of which being small local wavelengths (compared with the structures length) and unstable post-buckling responses.
In the present case, the shell wall in the vicinity of the crown appears to act as the rod, while the coupling of radial displacements and circumferential stretching provides the elastic foundation. As the end-shortening increases further, the amplitude of a single inward half-wave increased rapidly (the peak radial velocity is ~360 mm/s) and the analogy of the rod on an elastic foundation continues. In contrast to the hybrid behavior of constant stiffness cylinders, in which radial displacements are localized in the axial direction but distributed around the circumference [26], the variation of stiffness here results in localization in both directions.
The magnitude of the peak radial displacement at point B was -20.07 mm (w/h = 9.80), which is larger than the experimentally measured value of ~10 mm (w/h ≈ 5). This error is consistent with the differences in the predicted and measured load drops at this part of the loading cycle. Immediately after the peak load is reached, the radial displacement pattern changes from a non-localized harmonic one to a highly localized dimpled one. The peak radial velocity of the unstable inward half wave is ~430 mm/s. The peak radial displacement of the shape shown in Fig. 15c

Strain energy analysis
The force versus end-shortening response of Shell B differs from that of Shell A in two respects.
Firstly, a smaller drop in load was observed following the initial buckle of Shell B than Shell A, and secondly, the occurrence of the first localization in Shell B had a relatively small effect on its axial stiffness. Since both of these characteristics are interesting from a design point of view, a detailed comparison of their post-buckling behavior by means of a strain energy analysis was conducted. Figure 15. Deformed surface plots of Shell B at and immediately after the critical load is reached.
The following energies were computed from the element strains and stress-resultants using numerical integration (i.e. the rectangle rule): , , and , , and ; where: N x , N y and N xy are the membrane stress-resultants; M x , M y and M xy are the bending stressresultants; ε x , ε y and ε xy are engineering strains; κ x , κ y and κ xy are the changes in curvature and twist; and A is the total surface area of the shell. The total stored energy of the shells, neglecting energy stored due to transverse shear, is the sum of these components; that is:  are due to the existence of a boundary layer at each end of the shell. If the shell's ends were free to expand these components would be zero.
At the point of buckling, a portion of the prebuckling strain energy is converted into kinetic energy, resulting in a drop of the total strain energy. This is more pronounced in Shell A than in Shell B. In addition to the conversion into kinetic energy, the buckling event causes a redistribution of strain energy, from a single mode of storage, to one involving each of the stretching and bending components. In the case of Shell A, the strain energy stored due to membrane stretching reduces by approximately 68% and is redistributed into strain energy due to circumferential bending, axial bending and twisting (in order of decreasing magnitude). The propensity for the development of large quantities of strain energy due to circumferential bending is due to the formation of a series of 'knuckles' which separate inward dimples in the y-direction (see Figs. 11b & 14d). As the amount of end shortening increases, so does the energy stored due to bending, while the membrane energy appears to reach an upper limit.
In the case of Shell B, on the other hand, the formation of the first localization corresponds to a drop in axial stretching energy of only 33%. On inspection of the shells' deformed shapes immediately after buckling, it is clear that the first localization occurring in Shell B occupies a smaller part of the shell's circumference than that of Shell A. Coupling this information with the distributions Interestingly, one may draw parallels between such behavior and the principle of the effective width used by von Kármán [27] to explain the post-buckling behavior of plates. Here, a certain part of a plate's physical width is assumed unable to carry in-plane loads; due to the development of large transverse displacements. While its in-plane stiffness is degraded, the plate material in the proximity of its boundary continues to develop membrane stress, giving the structure post-buckling strength reserves. In the case of Shell B, the localization of radial displacements has resulted in a similar situation; albeit a more precarious one. The localization has resulted in a reduction in the shell's effective circumference, where the side regions carry the majority of the axial load. Unlike a plate boundary, however, the side regions are vulnerable to local buckling themselves, and collapse after continued loading.
The reasons for the differences in behavior are difficult to predict from the numerical analysis, however, there must be some dependency of the size of the localization on the overlaps in the crown/keel regions. In the case of Shell B the reduction of overlaps resulted in a laminate with smaller axial bending stiffness compared with Shell A. It seems likely, therefore, that for a relatively small and localized radial displacement field to occur, a sufficient local increase in membrane stiffness must be present (to attract compressive load and cause local buckling), but with a small enough local bending stiffness to result in the formation of short wavelengths.

Conclusions
The behavior of two variable stiffness cylindrical shells under axial compression has been investigated by linear and nonlinear finite-element analysis and compared with experimental results.
Details, such as the effects of the experimental boundary conditions, rotations of the test platens and variations in the ply thickness, were taken into account and discussed.
Linear and nonlinear predictions of the buckling loads were both found to give similar results which ranged in error from 0.4-7 % compared with the experiment. The effect of unintentional platen rotations which occurred in the experiment (between 10 and 20 kN/degree) was included in the linear analyses and was found to reduce the predicted buckling load by up to 2%.
Nonlinear analysis results were found to be in good agreement with the experimental measurements, in terms of global displacement patterns up to approximately 150% of the critical endshortening deformation. This accuracy made possible the visualization of the mechanism leading to collapse, which was found to be the rapid growth of an inward half-wavelength of the nonlinear buckling mode-shape.
A striking feature of the post-buckling behavior of these variable-stiffness shell structures was the presence of a stable and repeatable localization of the radial displacements, which forms after the initial buckling event. In the case of Shell A, the stiffness of the laminate was such that a localized dimple with a significant diameter compared to the cylinder's circumference formed, giving it a relatively low effective circumference and, hence, a low post-buckling stiffness. On the other hand, Shell B, which had a more uniform thickness, showed a first localization with a diameter which was relatively small in proportion to the shell's circumference, enabling a direct membrane load path through the structure. These differences in behavior were reflected in the components of the total strain energy which showed that the shell with a relatively large effective circumference stored the majority of its strain energy in membrane stretching.
These findings are interesting from a design perspective because they show that there is a sensitivity of a shell's post-buckling behavior on the precise form of its variable fiber-angle lay-up. In the future, structures with such effects could be exploited by tailoring the stiffness properties to give them novel, more benign, nonlinear post-buckling behavior.