Accepted Manuscript

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

Simon C. White, Paul M. Weaver, K. Chauncey Wu

PII: DOI:

Reference: To appear in:

S0263-8223(14)00669-2

http://dx.doi.org/10.1016/j.compstruct.2014.12.013 COST 6073

Composite Structures

Please cite this article as: White, S.C., Weaver, P.M., Chauncey Wu, K., Post-buckling analyses of variable-stiffness composite cylinders in axial compression, Composite Structures (2014), doi: http://dx.doi.org/10.1016Zi.compstruct. 2014.12.013

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

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

compression

Simon C. White"'1, Paul M. Weaver"'2 & K. Chauncey Wub 3

aAdvanced Composite Centre for Innovation and Science, University of Bristol, Queen's Building,

BS8 1TR

bStructural Mechanics and Concepts Branch, NASA Langley Research Center, Hampton VA

Abstract

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 postbuckled 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.

ficant am /

1 PhD Student. Email address: simon.white@bristol.ac.uk

2 Professor in Lightweight Structures. Email address: paul.weaver@bristol.ac.uk

3 Senior Aerospace Engineer, Structural Mechanics and Concepts Branch, RD. Email address: k.c.wu@nasa.gov

1. 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-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 Karman & 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 displacement-controlled 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, Huhne, 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.

2. 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 arc-length, y. We also use a Cartesian system (X,Y,Z) to define positions in three-dimensional space. Figure 1a shows the general arrangement of the specimens. Their nominal dimensions are length L = 889 mm and internal radius of curvature R = 207 mm. The shell walls are composite laminates with

Figure 1. (a) General arrangement of the shell (dimensions are in millimeters). (b) Potting and boundary hoop section detail.

nominal stacking sequences (defined in terms of the angles made between the x-axis and the fiber) designed to be of the form: [±45°,±^(y)]s. The function $(y) is a piecewise-continuous variable fiberangle which varies between the following design values:

/y = 0, -R ,-R,3-R J = 10° ,45° ,10° ,45°

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. D16 and D26) 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

Table 1. Specimen dimensions and tolerances

Shell A (overlaps) Shell B (no overlaps)

Mean Std. dev. Mean Std. dev.

Internal radius, R (mm) 206.86 0.53 206.91 0.48

Length, L (mm) 889.0 NA 889.0 NA

Laminate thickness, h (mm) 1.270 0.279 0.991 0.510

Average ply thickness, t (mm) 0.128 0.006 0.124 0.008

Shell A Shell B

(a) (b)

Figure 2. Courses of the variable angle plies generated by the AFP CAD software.

consequence, ply gaps and overlaps develop [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,103,-103]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/8552prepreg, 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.

Figure 3. Photo of the test rig.

In order to replicate a fully-clamped boundary condition, the ends of both shells were encased in epoxy, in between two aluminum hoops (see Fig. 1b). The ends of the assemblies were machined parallel and flat, resulting in an overall length of 889.0 mm and a free shell length of 838.2 mm. The axial compression force, designated X, was measured with a load cell, and axial displacements were measured with four displacement transducers, one attached to each corner of the test machine platens. 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).

3. 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.

3.1. 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 mm2 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.

>TM stand;

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 E11, E22, and G12 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).

3.2. Damping

The Hilber-Hughes-Taylor (HHT-a) [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-a 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 (a = -1/3) was required to reduce the number of time increments to a reasonable unt, or even in some cases achieve convergence of every time increment.

Figure 4. The finite-element mesh used in the analyses. The red regions represent those

contained within the encasement.

Table 2. Material properties.

Elastic moduli (GPa) Poisson's ratio Mass density (g/cm3)

Material Eu E22 G12 V12 P

IM7/8552 134.3 9.38 5.21 0.347 1.57

Epoxy 3.00 3.00 1.15 0.300 1.20

Aluminum 70.0 70.0 26.9 0.300 2.80

3.3. 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 constrained degrees of freedom were radial translations and rotations about the j-axis. This was represented in the model by restraining translations of the nodes in the encased region in the Y-Z plane. A further effect of the encasement was to provide an increase in the local axial stiffness of the shell wall which, in turn, increases the shell's global axial stiffness. While this may not have a large effect on the buckling loads, a poor estimate of the prebuckling stiffness affects the accuracy of the simulated load-deflection paths in a nonlinear analysis. This effect was investigated by Hilburger & 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.

Figure 5. (a) The finite-element mesh used in the boundary sub-model. (b) A plot of a typical deflected shape of the boundary sub-model.

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

(lX^ kenc (i.e., kenc

encased/kfree) was calculated as:

d\ d 2 dn

where d1 and d2 are nodal displacements (shown in Fig. 5b), and n 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: Ez = 9 GPa and Gxz = 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, n, 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.

Table 3. Encasement sub-model results for different stacking sequences.

Elastic moduli (GPa)

Poisson's ratios

Stiffness ratio

Axial stiffness increase (%)

Location E E Eyy Vxy Vxz n

Sides, both shells 18.36 18.36 0.758 0.117 2.88 3.87

Crown, Shell B 71.00 21.93 0.744 0.142 1.46 1.83

Crown, Shell A 97.32 16.51 0.728 0.159 1.16 0.81

The results show that the stiffness ratio, n, is inversely proportional to the laminate's axial elastic modulus and that its maximum value is at the sides of both shells. The increase in the shells' local axial stiffness due to the encasement is in the range of 1.46 < n < 2.88 for Shell A and 1.16 < n < 2.88 for Shell B. Inclusion of the encasement in general resulted in a small increase in the shells' axial stiffness (+3.87% maximum). This result is due to the relatively short axial length of the encased region compared to the shell's total length.

;ly short

. are pre

4. Results and discussion

Numerically predicted results for the two specimens are presented in this section. First, the results of a short investigation into the effects of variations in the encasement stiffness ratio and the ply thickness on the non-linear buckling predictions are presented. The effects of small platen rotations on the linear buckling predictions of both shells are then investigated. This was motivated by observations of platen rotations in the experiment. The remaining sections are concerned with the results of the non-linear collapse simulations. Full loading and unloading cycles are predicted and compared with the experimental data, the initial collapse mechanisms are explored in detail and, finally, differences in the behavior of the shells are investigated by means of an energy analysis.

4.1. 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 n 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 n equal to small (1), intermediate (2) and large (4) values (achieved by reducing the encased elements length by a factor 1/n). This relatively large range of n was chosen to produce a correspondingly large variation in the dependent variables (i.e. the pre-buckling 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 n was exhibited by the prebuckling stiffness (as much as 3.9% variation for shells A and B with n in the range 1 < n < 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 n =1) to 57.00% (with n = 2), a difference of 2.13%. Following a further end-shortening, however, the loads for all of the values of n became approximately equal. In the case of Shell B, these values were 27.27% (with n = 4) and 26.51% (with n = 1), a difference of 0.76%.

In general, the predictions of the behavior of Shell B were worse than those of Shell A. Errors in the prebuckling stiffness and buckling load (with n = 1) were -3.5 and -15.5 %, respectively. Given the relative magnitudes of these errors and noting that the axial stiffness and the buckling loads of an isotropic circular cylindrical shell scale with h and h2, respectively [4], it is likely that the discrepancies are related to an error in the value of the homogeneous ply thickness t*. Indeed, since the buckling load is proportional to sixty-four times the square of the individual ply thickness one can expect a large error to result from a relatively small error in ply thickness.

Figure 7a shows measured frequency of ply thickness for Shell B which has been calculated by dividing the laminate measurements by eight. There is a small but clear skewing in the distribution

Figure 6. Initial buckling results for: low (1), intermediate (2) and high (4) values of n. (a) Results for Shell A (t* = 0.128 mm). (b) Results for Shell B (t* = 0.124 mm).

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-buc 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% in crease 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 the intermediate stiffness ratio n = 2 has been used for both shells.

s of t*

t* = 0. 128

. 128 mm and

4.2. 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 and 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.

Figure 8 shows the rotation measurements from the experiments. In the case of Shell A, the forcerotation responses about both the Y- and Z-axes were fitted relatively well with line are regressions

Figure 7. (a) Measured ply thickness frequency (1/8th the measured laminate thickness). (b) Load vs. end-shortening results for Shell B with t* = 0.128 mm.

having slopes of approximately 18.98 (R2 = 0.987) and 14.04 (R2 = 0.992) kN/degree. These are shown by dashed lines in Fig. 8a. In the case of Shell B, on the other hand, only rotation about the Z-axis showed an acceptable linear force vs. rotation trend. A least squares fit gave a slope of 9.73 (R2 = 0.719) kN/degree. Rotation about the 7-axis, in this case, showed an initially linear response but only over a small part of the prebuckling step.

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.

Table 4 shows the buckling eigenvalues obtained with various combinations of platen rotations included in the load case. Note that in the experiments the buckling loads were measured as 172.64 and 76.46 kN [18]. Predictions of the buckling loads where platen rotations were excluded have an error compared with the experimental values of -0.38 and -8.19% for Shells A and B, respectively. These errors are smaller than those reported in ref. [ 18] because in the present work the ply material

Figure 8. Measured relative rotations of the test machine platens. (a) Results for Shell

A. (b) Results for Shell B.

Table 4. Buckling eigenvalue results (kN) for Shell A with various combinations of applied

platen rotation.

Applied rotations

Model Zero rotation ^Z and

Shell A, perfect 172.42 168.98 172.09 169.06

Shell A, imperfect 171.99 168.12 171.63 169.03

Shell B, perfect 71.647 69.59 NA NA

Shell B, imperfect 70.20 68.48 NA NA

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).

Figure 9 shows the linear buckling modes for the two shells with the applied platen rotations. In contrast to the modes found without platen rotations (see ref. [18]), the buckling modes are concentrated toward the crown and have larger radial displacements near the end X = L. This effect is due to the positive bending moments about the Z axis increasing the compressive axial membrane stress in the crown. The bucking mode of Shell A is wider in the circumferential direction than that of Shell B and has a smaller axial wave number in the x-direction (14 compared with 20). This behavior is due to the larger bending stiffness of the crown laminates in Shell A than Shell B in both the x- and j-directions.

liscussed. A

es in Shel

4.3. 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

Figure 9. Linear buckling modes under combined end-shortening and end-rotation. (a) Results for Shell A with imperfect geometry. (b) Results for Shell B with imperfect

geometry.

(see Fig. 8). Platen rotations about the Z- and 7-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.

4.3.1. Shell A

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].

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 shows a comparison of the experimental and simulated radial displacements at the initial point of buckling (point A, Fig. 10b). The most significant displacements were confined to the crown/keel regions where the highest compressive stresses were found. The number of axial half-waves in the experiment and the simulation are 19 and 17, respectively.

End shortening beyond the point at which the peak load developed resulted in the formation of an

Figure 10. Load vs. end-shortening results for Shell A based on a model with perfect (a) and imperfect (b) cylindrical mesh geometries.

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

ftpp.hr i: ,

Ahaqus

(a) Point A, buckling

(b) Point B. post-buckling

(c) Point C, deep post-buckling

(d) Point D, deep post-buckling

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.

4.3.2. Shell B

Figure 12 shows a comparison of the force versus end-shortening results obtained fri experiment and the simulation. Figure 13 shows qualitative color maps of radial dis (experiment) and displacement magnitudes (simulation), plotted on the deformed st different points on the loading cycle.

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

Figure 12. Load vs. end-shortening results for Shell B based on a model with perfect (a) and imperfect (b) cylindrical mesh geometries.

section 4.5). The predicted shell axial stiffness after the formation of the first localization is in good agreement with the experimental observations (-2.12 % error).

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

(a) Point A, buckling

(b) Point B. post-buckling

(c) Point C, deep post-buckling

(d) Point D, deep post-buckling

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.

4.4. 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 ^'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.

e 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.

Figure 15 shows the incremental changes in Shell B's deformed geometry, again, occurring immediately after the buckling load (i.e. at points on the line A-B of Fig. 12b). At the point of buckling (point A, Fig. 12b), a harmonic radial displacement pattern focused at the crown forms with a maximum radial displacement of 0.41 mm (w/h = 0.39). These deformations are larger than the experimentally measured values of ~0.10 mm (w/h = 0.1), yet with the same order of magnitude. 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 is 10.04 mm (w/h = 9.80). After a small time period the central dimple stabilizes without an appreciable increase in end shortening (see Fig. 13b) with an amplitude of 15.24 mm (w/h = 14.88). Again, localization in both the axial and circumferential directions is observed. The peak experimental radial displacement at this stage is 10.16 mm (w/h = 9.92). This difference is consistent with the fact that the predicted load drop is larger than the experimental value.

The buckling behavior of these structures differs significantly from that of traditional constant stiffness structures. In the post-buckling range, we also see interesting new force-displacement behavior, particularly in the case of Shell B. In the next section their relative post-buckling behavior is explored further by means of a comparative strain energy analysis.

Figure 15. Deformed surface plots of Shell B at and immediately after the critical load

is reached.

4.5. 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.

The following energies were computed from the element strains and stress-resultants using numerical integration (i.e. the rectangle rule):

- J NxexdA, = - J NysydA, ^ = - J Ny^

where: Nx, Ny and Nxy are the membrane stress-resultants; Mx, My and Mxy are ti

Bx = - J MxKxdA, By = - J My ^ydA, and By = - J i^^dA ;

stress-and twist; eglecting energy

resultants; sx, sy and sxy are engineering strains; kx, ky and kxy are the changes i and A is the total surface area of the shell. The total stored energy of the si stored due to transverse shear, is the sum of these components; that is:

E = S + Sy + Sxy )+(Bx + By + Bxy )

Figures 16 and 17 show the time-wise development of stored strain energy in the shells, plotted in terms of the individual membrane and bending energy components. In the prebuckling regime both structures store the majority of their strain energy due to axial compression, with the exception of small quantities stored due to axial bending and circumferential stretching. These latter components

Figure 16. Strain energy components of Shell A over the course of the loading step.

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 j-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 Coupling this information with the distributions

120r Pre-buckling

Post-buckling

0 100 200 300 400 500 600

Time (s)

Figure 17. Strain energy components of Shell B over the course of the loading step.

of energy components illustrated in Figs. 16 and 17, it is evident that Shell B has buckled into a configuration which maintains a load path along the shell's sides and that it is based on the development of membrane stress. This explains the retention of axial stiffness into the post-buckling range.

Interestingly, one may draw parallels between such behavior and the principle of the effective width used by von Karman [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

ormation

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 end-shortening 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

et Stanford

6. Acknowledgements

The authors would like to thank Dr. Bret Stanford for his help in the generation of the finite-element models. The authors would also like to thank the Engineering and Physical Sciences Research Council for supporting the ACCIS Centre for Doctoral Training. Grant No. EP/G036772/1.

7. References

[1] Hühne, C., Rolfes, R., Breitbach, E., & Teßmer, J. (2008). Robust design of composite cylindrical shells under axial compression—simulation and validation. Thin-Walled Structures, 46(7), 947-962.

[2] Degenhardt, R., Kling, A., Bethge, A., Orf, J., Kärger, L., Zimmermann, R., & Calvi, A. (2010). Investigations on imperfection sensitivity and deduction of improved knock-down factors for unstiffened CFRP cylindrical shells. Composite Structures, 92(8), 1939-1946.

Bisagni, C., & Cordisco, P. (2003). An experimental investigation into the buckling and post-buckling of CFRP shells under combined axial and torsion loading. Composite Structures, 60(4), 391-402.

[4] Calladine, C. R. (1989). Theory of shell structures. Cambridge University Press.

[5] Gerard, G. & Becker, H., (1957). Handbook of structural stability part III - The buckling of curved plates and shells. N.A.C.A. T.N. 3783.

[6] Simitses, G. J. (1986). Buckling and postbuckling of imperfect cylindrical shells: a review. Applied Mechanics Reviews, 39(10), 1517-1524.

[7] Almroth, B. O. (1965) Influence of edge conditions on the stability of axially compressed cylindrical shells. N.A.S.A. CR-161.

[8] von Kärmän, T., & Tsien, H. S. (1941). The buckling of thin cylindrical shells under axial compression. Journal of Spacecraft and Rockets, 40(6), 898-907.

[9] Koiter, W. T., On the Stability of Elastic Equilibrium. (in Dutch), H. J. Paris, Amsterdam, Holland, 1945; translation available as AFFDL-TR-70-25, February, 1970, Wright-Patterson Air Force Base.

[10] Arbocz, J. (2000). The effect of imperfect boundary conditions on the collapse behavior of anisotropic shells. International Journal of Solids and Structures, 37(46), 6891-6915.

[11] Hilburger, M. W., & Starnes Jr, J. H. (2002). Effects of imperfections on the buckling response of compression-loaded composite shells. International Journal of Non-Linear Mechanics, 37(4), 623-643.

[12] Hilburger, M. W., & Starnes Jr, J. H. (2004). Effects of imperfections response of composite shells. Thin-Walled Structures, 42(3), 369-397.

[13] Lancaster, E. R., Calladine, C. R., & Palmer, S. C. (2000). Paradoxical buck a thin cylindrical shell under axial compression. International Journal Sciences, 42(5), 843-865.

ms of the f

ickling beh

uckling

behaviour of of Mechanical

Guggenberger, W., Greiner, R., & Rotter, J. M. (2000). The behaviour of locally-supported cylindrical shells: unstiffened shells. Journal of Constructional Steel Research, 56(2), 175-197.

Zhu, E., Mandal, P., & Calladine, C. R. (2002). Buckling of thin cylindrical shells: an attempt to resolve a paradox. International Journal of Mechanical Sciences, 44(8), 1583-1601.

Calladine, C. R. (2002). A shell-buckling paradox resolved. In Advances in the Mechanics of Plates and Shells (pp. 119-134). Springer Netherlands.

Deml, M., & Wunderlich, W. (1997). Direct evaluation of the 'worst' imperfection shape in shell buckling. Computer Methods in Applied Mechanics and Engineering, 149(1), 201-222.

Wu, K. C., Stanford, B. K., Hrinda, G. A., Wang, Z., Martin, R. A., & Kim, H. A. (2013). Structural Assessment of Advanced Composite Tow-Steered Shells. In Proceedings of the 54th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference. Boston, Massachusetts. Paper AIAA-2013-1769.

Wu, K. C., Tatting, B. F., Smith, B. H., Stevens, R. S., Occhipinti, G. P., Swift, J. B., & Thornburgh, R. P. (2009). Design and Manufacturing of Tow-Steered Composite Shells Using Fiber Placement. In Proceedings of the 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference. Palm Springs, California. Paper AIAA-2009-2700.

[20] Wu, K. C. (2008). Design and analysis of tow-steered composite shells using fiber placement. In Proceedings of the American Society for Composites 23rd Technical Conference. Memphis, TN.

[21] Lukaszewicz, D. H. J., Ward, C., & Potter, K. D. (2012). The engineering aspects of automated prepreg layup: History, present and future. Composites Part B: Engineering, 43(3), 997-1009.

HexPly® IM7/8552 Product Data Sheet, Feb 2013. Hexcel Corp.

Dassault Systèmes (2012) Abaqus/Standard (version 6.12-1). URL: www.3ds.com/products-services/simulia/portfolio/abaqus/latest-release/

Hilber, H. M., T. J. R. Hughes, and R. L. Taylor (1977). Improved Numerical Dissipation for Time Integration Algorithms in Structural Dynamics. Earthquake Engineering and Structural Dynamics, vol. 5, pp. 283-292.

Hunt, G. W., Bolt, H. M., & Thompson, J. M. T. (1989). Structural localization phenomena and the dynamical phase-space analogy. Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 425(1869), 245-267.

[26] Hunt, G.W., Lord, G.J., & Peletier, M.A. (2002). Cylindrical shell buckling: a characterization of localization and periodicity. Discrete and Continuous Dynamical Systems - Series B, vol. 3(4), 505-518.

[27] von Karman, T., Sechler, E. E., & Donnell, L. H. (1932). The strength of thin plates in compression. Trans. ASME, 54(2), 53-57.