## Abstract

This study introduces an approximate analytical model to predict the post-buckling response of cylinders with tailored non-uniform distributed stiffness. The shell's wall thickness, and thus its stiffness, is tailored so as to obtain multiple controlled elastic local buckling events when the cylinder is subjected to uniform axial compression. The proposed model treats cylinder segments of different stiffness as individual panels and combines their response by considering them as connected linear or nonlinear springs. The governing equations for the panels are formulated using von Karman's theory and solved by Galerkin's approximate method for a predefined radial deformation. Radial deformation functions are used to improve the model's accuracy, and results show that the model's accuracy increases significantly with the number of considered radial functions. The model's predicted axial response for different cylinders is compared with results from experiments on three-dimensional printed samples. Results indicate that this model accurately predicts the order of the buckling events, while the buckling forces from the model are higher than those measured experimentally.

## 1 Introduction

Thin-walled cylindrical shells are widely used in many engineering fields, such as offshore structures and aerospace vehicles due to their high load-carrying capacity despite their lightweight. Yet, their slender wall thickness makes them highly susceptible to buckling when subjected to axially compressive loads. Engineers have studied the buckling strength of cylindrical shells for many years in order to predict and thus avoid such behavior. However, because of the high sensitivity of buckling strength to initial imperfections, the discussion on determining design factors for the axial strength of cylindrical shells is still open [1]. On the other hand, the elastic post-buckling response of structures has garnered significant research interest in recent years due to newly emerged smart applications arising from such behavior for energy dissipation, energy harvesting, smart switches, and self-powered sensors [210]. However, predicting and controlling the elastic post-buckling response of thin-walled cylindrical shells is even more challenging than predicting their buckling strength. It has been shown that if initial imperfections are simulated precisely, the finite element (FE) method is able to accurately predict the buckling and post-buckling response of cylindrical shells [11]. Since it is not possible to determine the initial imperfections of a cylinder prior to design and manufacturing, researchers continue to look for reliable methods to predict the buckling response in a way that considers the most probable initial imperfections [1215]. However, despite the numerous methods to consider initial imperfections in analytical and numerical models, none of them can guarantee the most efficient and reliable design of cylindrical shells [16].

More recently, researchers have tried to take control over the buckling and post-buckling response of axially loaded cylindrical shells by manipulating the cylinder's shape to trigger specific buckling modes. Ning and Pellegrino [17] showed that changing the cross section of the cylinder by applying smooth waves along its circumference reduces uncertainty on their buckling strength in addition to increase the buckling critical load. Cox et al. [18] used seeded modal imperfections to control the post-buckling response of the cylindrical panels. Using finite element analyses (FEA), they showed that such seeded imperfection has a significant positive effect on controlling the post-buckling response of the cylindrical panels without a significant increase in the mass of the shell. In a series of studies, Hu and Burgueño proposed different types of seeded imperfections on cylindrical shells to control their post-buckling response [1923]. The approaches included changing the pattern of the shells’ wall stiffness, the provision of external and internal radial constraints, and the seeding of sinusoidal imperfections to the cylinder's wall. Their results showed that seeded modal imperfections or providing a patterned non-uniform wall stiffness distribution significantly decreased the sensitivity of the post-buckling response to initial imperfections. The studies by Hu and Burgueño further showed that providing a tailored non-uniform stiffness distribution (NSD) to the cylinder's wall (NSD cylinder), achieved by thickening of different regions of the cylinder (see Fig. 1), allowed control of the number, location, and sequence of local buckling events in the cylindrical shell [23]. In follow-up studies for this type of cylinder design, Guo et al. [24] presented a model that described the effect of localized wall thickening on the strain energy released due to the localized buckling events, which is manifested as load drops in the cylinder's global force-deformation response. However, the approach did not permit predicting the exact force, displacement, and force drop for a given thickening pattern.

Fig. 1
Fig. 1
Close modal

In this paper, a semi-analytical model that is able to predict the post-buckling response of axially loaded NSD cylinders is presented. The motivation to develop this model includes the following: (1) to improve an understanding of the parameters controlling the response and design of NSD cylinders, (2) to avoid the meshing, element, and discretization considerations in the finite element method, and (3) to facilitate design optimization. The model follows three steps for any given NSD cylinder: (1) dividing cylinder to parallel axially compressed segments, (2) predicting the response of each cylindrical panel, and (3) integrating the response from the individual panels. The stiffness is non-uniformly distributed on the wall of each segment of the NSD cylinder. Some panels are stiffer at both ends in the axial direction, and others are stiffer in the center.

The problem of instability in cylindrical shells with variable thickness has been studied by several researchers. Koiter et al. [25] investigated the buckling of a cylinder with harmonic variation in its thickness along the axial and circumferential directions of the cylinder, which is not the case in our study. Chen et al. [26] investigated cylindrical shells with stepwise thickness distribution under radial pressure. They kept the length of the cylinder constant and defined an equivalent thickness in a way that the stored energy after buckling in the equivalent cylinder matched the strain energy in the original cylinder. Using this method, they did not match axial stiffness of the cylinder since it does not play a key role in the radial pressure situation. Several studies have developed analytical solutions for axially compressed cylindrical shells with small stepwise variations in thickness only in the axial direction [2729]. In these studies, both thick and thin parts of the cylinder experience nonlinear deformations. However, in the current study, the difference in wall thickness of thick and thin regions is larger than 100% of the smallest thickness. As a result, buckling occurs only in the thin-walled areas. Also, the wall thickness of the panels varies in both the axial and circumferential directions.

For a cylinder subjected to radial pressure with stepwise thickness variation, Greiner [30] assigned an equivalent length for the part of the cylinder that experienced buckling. The equivalent length depended on the stiffness of the remaining part of the cylinder. In this work, the equivalent length approach is adopted to replace the thicker zones at the top and bottom of a panel by shell regions of equal thickness than that at the panel's mid-section. The equivalent length is calculated in a way that the lateral stiffness at the boundary of the thin areas does not change. The post-buckling response of an individual panel is formulated based on classic shell theory and solved with Galerkin's approximation. The approach can be used for any pattern of thickened areas; but in this work, it is used for the thickened areas at that separate the potential buckling panels, namely, top, bottom, and middle regions along the shell (see Fig. 1).

Section 2 discusses the division of the cylinder into parallel segments and the determination of the equivalent length of cylindrical panels with variable thickness. In Sec. 3, a semi-analytical predictive model for post-buckling response of an axially compressed cylindrical panel with clamped boundary conditions is presented. Validation of the proposed model is also presented in this section. Using the proposed model in Sec. 3 and the simplifications in Sec. 2, the post-buckling response of four NSD cylinders is predicted and compared with the response from experiments in Sec. 4. Section 4.2 presents a parametric study on post-buckling response of axially compressed cylindrical panels conducted for panels with varying dimensions, and design contours are provided for NSD cylinders with a specific radius and body thickness.

## 2 Modeling Concept and Hypothesis

In NSD cylinders, buckling events are limited to slender zones (areas with thickness equal to tb in Fig. 1), where thicker zones (areas with thickness equal to t2 in Fig. 1) impede interaction between the slender zones. It is assumed that the parallel zones of NSD cylinders are isolated from each other and thus the cylinder is divided into parallel panels with clamped boundary conditions (C) along all sides. Figure 2 shows the separated parts of the cylinder and their boundary conditions. The side boundaries of panels that are subjected to buckling (panels A, B, and C in Fig. 2) are connected to the thick areas in panels D.

Fig. 2
Fig. 2
Close modal

Here, it is assumed that side boundaries are clamped but can move along the axial direction. Panels D with a very small thin shell region are not expected to experience buckling and as a result can be considered as linear axial springs. The separated panels (A, B, C, and three panels D) are parallel to the loading direction, and since a uniform deformation is applied to the top edge of the cylinder, all these panels experience the same deformations at their top. Thus, the total axial load in the cylinder is the sum of forces in the panels, and they can be modeled as parallel springs.

Panels with buckling potential regions (panels A, B, and C in Fig. 2) have two thicker zones at the top and bottom and a thin zone at the center, which is the targeted buckling region. Thus, the main problem is to analyze the post-buckling response of a cylindrical panel with non-uniform thickness under compression. To simplify the panel with non-uniform thickness, the wall thickness of the top and bottom zones is made equal to the wall thickness of the center zone by modifying its length. The equivalent length of the thick-walled region at the top and bottom of the panel should be determined in a way that lateral stiffness at the edge of the thin areas is unchanged. The equivalent length of the shell with transformed thickness (top and bottom regions) is calculated in two general steps. First, the thicker regions at the top of the panels are assumed as one-way retrained cantilevers and their tip transverse displacement, Δ, under an applied moment, M, is determined (see Fig. 3). Second, determine the length of a cantilever shell with the thickness of the mid-section of the panel (tb in Fig. 3) that yields the same deformation Δ under equal bending moment M.

Fig. 3
Fig. 3
Close modal
The bending moment and deformation at the edge of the cantilever elements are related by Eq. (1):
$M=3EI1il1i2Δ=3EIbl1i′2Δ$
(1)
where i indicates the panel's name (A, B, or C), E is the modulus of elasticity, $l1i$ is the initial length of thicker areas at the top or bottom, $l1i′$ is the equivalent length of the same area, and tb and $t1i$ are the thicknesses of the panel's mid-section and the thicker areas, M is the bending moment of the boundary, and Δ is the deformation of the edge in the radial direction. The equivalent length of the thicker zones is then calculated with Eq. (2):
$I1il1i=I1i′l1i′⇒l1i′=(tbt1i)3/2l1i$
(2)
The cylindrical panels can now be modeled with a uniform thickness, which simplifies determining an analytical solution to the post-buckling response of each panel. However, the equivalent shell panel has a different axial stiffness than the original one. To correct the axial stiffness deviation, a linear spring is added in series to the equivalent panel so as to match initial panel stiffness. The stiffness of this spring is calculated with Eq. (3):
$ki1=bEt1itb(b1t1i+btb−b1tb)btb((l1i′−l2)t1i−l1itb)−b1(t1i−tb)(l1i′t1i+l1itb)$
(3)

The total axial displacement in an individual panel is equal to the sum of the displacements from the linear springs and the displacement from the nonlinear response of the panel. Similarly, the axial force corresponding to a given displacement in the cylinder is the sum of axial forces in all panels (A, B, C, and three D panels) at that displacement. In the other words, the panels contribute to the response of the cylinder as resisting elements in parallel, while different regions (i.e., different thickness) within an individual panel behave as elements in series. Figure 4 shows the framework followed in this paper to find the analytical solution for the post-buckling response of the NSD cylindrical shell in Fig. 1.

Fig. 4
Fig. 4
Close modal

The stiffness of all the equivalent linear springs in the cylinder are known, except the axial nonlinear response of the equivalent cylindrical panels $KAnl$, $KBnl$, and $KCnl$ in Fig. 4. In Sec. 3, a classical method is used to develop a relation for the nonlinear response of an individual cylindrical panel with all-clamped boundary conditions.

## 3 Semi-Analytical Model for Axially Compressed Cylindrical Panels

### 3.1 Analytical Solution for Cylindrical Panel With All Edges Clamped.

In Sec. 2, the concepts of separating panels and transforming them into equivalent panels with uniform thickness were discussed. This section presents the analytical solution for the post-buckling response of a single panel with clamped boundaries at the top and bottom and a guided roller on the sides. The resulting force–displacement curves are then compared with the results from FE analyses for individual panels with equal boundary conditions. Figure 5(a) shows the panel, the boundary conditions, and the direction of the applied compressive displacement. Notations for the coordinate system, deformations, and stresses in each direction for the analytical and FE models are shown in Fig. 5(b).

Fig. 5
Fig. 5
Close modal
According to classical shell theory, the deformation of cylindrical shells should satisfy the following equilibrium and compatibility equations based on von Karman's theory [31]:
$D∇4w−f,xxR−[f,yy(w,xx+w,xx*)−2f,xy(w,xy+w,xy*)+f,xx(w,yy+w,yy*)]=0(∇4=∂4∂4x+2∂4∂2x∂2y+∂4∂4y),(,x=∂∂x,,y=∂∂y)$
(4)
$∇4f−Et(w,xy2−w,xxw,yy−w,xxR+2w,xyw,yy*−w,xxw,yy*−w,yyw,xx*)=0$
(5)
In these equations, R is the radius of the panel, w is the radial deformation, w* is the initial radial imperfection, and f is the airy stress function. Based on the definition of the airy stress function, the stress in the x-direction, stress in y-direction, and shear in the xy surface are given by
$Nx=f,yy,Ny=f,xx,Nxy=−f,xy$
(6)
The in-plane deformations u and v are related to out-of-plane deformation w and w* and the airy stress function if the following way:
$f,yy=Et1−ν2[u,x+12w,x2+w,xw,x*+ν(v,y−wR+12w,y2+w,yw,y*)]f,xx=Et1−ν2[ν(u,x+12w,x2+w,xw,x*)+v,y−wR+12w,y2+w,yw,y*]−f,xy=Et2(1+ν)[u,y+v,x+w,xw,y+w,x*w,y+w,xw,y*]$
(7)
The solution for Eqs. (4) and (5) is based on taking the radial deformation w and the radial initial imperfection w* as functions of coordinates x and y with unknown coefficients in a way that they satisfy the boundary conditions of the panel. It is then possible to find the unknown coefficients by minimizing the integration of Eqs. (4) and (5) with respect to them. Considering the clamped boundary condition in all edges with the panel sides allowed to move axially, w and w* should satisfy Eq. (8):
$w=∂w∂x=w*=∂w*∂x=0(forx=0andl)w=w*=0(fory=0andb)fxy=0(fory=0andb)$
(8)
Following radial displacement and initial imperfections can satisfy the boundary conditions:
$w=∑i1=1Λ∑j1=1ψ∑i2=1Λ∑j2=1ψai1i2j1j2(sini1πxl)(sini2πxl)(sinj1πyb)×(sinj2πyb)w*=∑i1=1Λ∑j1=1ψ∑i2=1Λ∑j2=1ψμi1i2j1j2t(sini1πxl)(sini2πxl)(sinj1πyb)×(sinj2πyb)$
(9)

Here, x and y are the in-plane coordinates of the panel, i1, i2, j1, and j2 are the number of the waves in the x- and y-directions (respectively), ai1i2j1j2 is the unknown coefficient for the radial deformation, μi1i2j1j2 is the known coefficient for the initial imperfection, and t is the panel thickness. A larger number of modes (Λ) considered in the solution is expected to provide a more accurate force–displacement response. However, increasing the number of modes in the deformation functions increases the computational cost of the solution.

By replacing w and w* from Eq. (9), Eq. (5) can be written as follows:
$∇4f=∑i=0Λ(Λ+1)∑j=0ψ(ψ+1)∑k=12fijkcos(iπxl+j(−1)kπyb)$
(10)
where Λ and ψ are the maximum number of contributed modes in Eq. (9) for the x- and y-directions, and fijk is the multiplier for the cosine function in ∇4f. Solving Eq. (10) results in a function F that has similar cosine functions to ∇4f but different multipliers to those cosine expressions (Fijk)
$F(x,y)=∑i=0Λ(Λ+1)∑j=0ψ(ψ+1)∑k=12Fijkcos(iπxl+j(−1)kπyb)$
(11)
However, the airy stress function can also contain functions λ(y) and θ(x) with ∇4λ(y) = ∇4θ(x) = 0. Then, it can be expressed as follows:
$f=λ(y)+θ(x)+F(x,y)$
(12)
Since the side edges of the panel are movable in the direction of the applied axial force (see Fig. 5(a)) the integration of Nx over width of the panel should be equal to the applied force (P). Thus,
$−P=∫0bf,yydy⇒λ(y)=−(12)(Pb)y2$
(13)
Also, the moveable edges only can move parallel to the x-direction, and as a result, the change in distance between these edges should be zero during the loading process. This condition is fulfilled in the average sense as [32]
$V=∫0l∫0bv,ydydx=0$
(14)
where
$v,y=Ny−νNxEt+wR−12w,y2−w,y*w,y*$
(15)
Thus, θ(x) can be expressed as
$θ(x)==−ν(12)(Pb)x2$
(16)
Having the airy stress function related to the radial deformation allows solving Eq. (4) numerically for any axial deformation. The axial deformation itself is calculated from Eq. (17):
$δ=∫0lu,xdxu,x=Nx−νNyEt−12w,x2−w,x*w,x*$
(17)

By replacing f and w in Eq. (4) and using Galerkin's approximate method for each unknown in Eq. (9), a nonlinear equation is obtained. If Λ = ψ = 1, there is only one equation to solve, and if Λ = ψ = 4, there are 100 nonlinear equations, which is equal to the number of unknowns in Eq. (9). These equations are solved using Newton's method for small iterations of the mean square deflection of the cylindrical panel (integration of the square of the radial deformation over the panel surface divided by the panel area), and the ai1i2j1j2 coefficients are found. The mean square deflection, unlike the maximum radial deformation, constantly increases during the pre- and post-buckling response, and therefore, the numerical solution converges well at bifurcation points. This method can be used for different boundary conditions by modifying the equation for radial deformation.

### 3.2 Model Verification for Post-Buckling Response.

Four finite element models with different geometries were prepared using the software abaqus [33] to verify the semi-analytical model. The geometry of the panels is shown in Table 1. Four node shell elements with reduced integration points (S4R) were used to represent the cylindrical panels, and a linear elastic material was assigned with an elastic modulus of 1300 N/mm2 and Poisson's ratio of 0.2. The top and bottom edges of the panels are clamped, while the side boundaries are fixed for radial and circumferential deformations and outward rotations but deform in the axial direction. Loading was in the form of a uniform axial displacement, equal to 1.5% of panel's length, applied to one end of the panel. To eliminate the effect of vibrations from the buckling events, a dynamic implicit solver with quasi-static application was used with a maximum time-step of 0.01 s.

Table 1

Dimensions of cylindrical panels for the verification of the analytical model

mmmmmmmm
(a)8080500.5
(b)8040500.5
(c)4060500.5
(d)6040500.5
mmmmmmmm
(a)8080500.5
(b)8040500.5
(c)4060500.5
(d)6040500.5

There are two options for applying simulated initial imperfections: (1) conducting an eigenvalue linear analysis and importing modes of interest to the post-buckling model or (2) considering the imperfection at the node's coordinates. Here, initial imperfections were seeded by perturbing the model's nodal coordinates based on the first buckling mode with μ1111 = 0.05 and all other μi1i2j1j2 in Eq. (9) equal to zero. The panel edges were fully clamped with the ability of axial motion along y = 0, y = b, and x = l (see Fig. 5). A uniform axial compressive displacement equal to 0.015h was applied to the top edge of the panel.

For each panel in Table 1, four semi-analytical models were developed for Λ = ψ = 1 to Λ = ψ = 4 in Eq. (9). The material properties and initial imperfections in the semi-analytical models were identical to those in the finite element models.

Figure 6 shows that analytical model's accuracy improves as the number of modes contributing to the radial deformation equation increases, especially for larger cylindrical panels. That is, more modes should be used in the analytical model when evaluating bigger panels in order to obtain an accurate prediction of the post-buckling response. For all cases, the analytical model with Λ = ψ = 4 accurately predicts the displacement corresponding to the first buckling event. However, the panel's post-buckling stiffness predicted by the analytical model is higher than the response calculated by the FE model.

Fig. 6
Fig. 6
Close modal

While increasing the number of involved modes in the radial deformation (Λ and ψ) leads to improved accuracy, it also grows the analysis time since the number of nonlinear equations escalates significantly. For example, by increasing Λ and ψ from 4 to 5, the number of nonlinear equations to solve rises from 100 to 225, which significantly increases computational time.

## 4 Elastic Post-Buckling Response of Axially Compressed Non-Uniform Stiffness Cylinders

### 4.1 Verification of Modeling Framework.

To verify the proposed predictive model for post-buckling response presented in Secs. 2 and 3, its predictive results are compared with test data from four prototype non-uniform stiffness (NSD) cylinders under axial compression. The geometry of the cylinders, with reference to Fig. 1, is presented in Table 2.

Table 2

Geometry of NSD cylinders tested under uniform axial compression

Test unittb (mm)t1A (mm)t1B (mm)t1C (mm)t2 (mm)l1 (mm)$l1A′$ (mm)$l1B′$ (mm)$l1C′$ (mm)l3 (mm)b1 (mm)b2 (mm)
CYL-10.511.3413010.567.151.33603565
CYL-20.511.5413010.565.751.33603565
CYL-30.512413010.563.751.33603565
CYL-40.512.5413010.562.651.33603565
Test unittb (mm)t1A (mm)t1B (mm)t1C (mm)t2 (mm)l1 (mm)$l1A′$ (mm)$l1B′$ (mm)$l1C′$ (mm)l3 (mm)b1 (mm)b2 (mm)
CYL-10.511.3413010.567.151.33603565
CYL-20.511.5413010.565.751.33603565
CYL-30.512413010.563.751.33603565
CYL-40.512.5413010.562.651.33603565

The NSD cylinders were three-dimensional (3D) printed from acrylonitrile butadiene styrene (ABS+) material using a Fortus 250MC printer (Stratasys Ltd., Eden Prairie, MN), see Fig. 7(a). All cylinders had an effective length of 100 mm and an inner radius of 50 mm. The cylinders were provided with 10 mm wide thickened edges at the top and bottom for connection to test platens to approximate clamped boundary conditions (the platens have grooves to fit the cylinder's edges and constrain rotation and radial displacements). Also, to avoid stress concentrations at the interfaces between thicker and thinner areas in the main cylinder body, the thickened areas had a stepwise change in thickness. The tests were conducted on a universal testing frame (Instron 5982) under displacement control as shown in Fig. 7(b). Loading consisted of a uniform quasi-static ramp of axial shortening up to 0.5 mm.

Fig. 7
Fig. 7
Close modal

For each of the cylinders in Table 2, a finite element model was prepared using the software abaqus. The modeling details were as described in Sec. 3.2, with two differences. First, the initial imperfection was simulated by using the cylinder's first buckling mode with an amplitude of 0.05t. Second, the elastic modulus for the analytical model was determined by matching the initial stiffness of the experiment and finite element models, which resulted in elastic modulus equal to 1300 N/mm2.

The force-deformation post-buckling response for each cylinder in Table 2 was predicted by using the proposed semi-analytical model presented in Secs. 2 and 3 considering Λ and ψ equal to 4. The framework presented in Fig. 4 was followed, whereby the NSD cylinder is separated into panels and the equivalent length of each panel is calculated from Eq. (2), and post-buckling response of the A, B, and C equivalent panels is calculated as presented in Sec. 3. The response of the parallel system of linear and nonlinear is determined by adding the forces from all springs for equal displacement increments to the system. The axial post-buckling responses of the cylinders from experiments, finite element models, and the proposed semi-analytical model are plotted in Fig. 8.

Fig. 8
Fig. 8
Close modal

Comparison of the responses in Fig. 8 shows that the difference between corresponding displacements of the buckling events between the analytical model, FE model, and experiment is smaller than 0.02 mm, which is less than 4% of the total deformation. Interestingly, in general, the analytical model better matches the experimental data for the displacement at buckling events. For all cases, the response from the analytical model is stiffer than the experimental response and finite element solution. This is consistent with the results presented in Fig. 6 that also show the analytical post-buckling response of cylindrical panels being stiffer than that obtained from the FE model.

It can be observed that the experimental traces in Fig. 8 initially show a lower stiffness with a stiffening response after approximately 0.1 mm of total deformation. This behavior is primarily attributed to test setup settlement, that is, the imperfect contact between the contact surfaces during loading. Another source for this response can be manufacturing imperfections, particularly at the transitions between thin and thick regions in the cylinder wall.

It can also be observed that the finite element analyses predict the magnitude of the force drops at the buckling events much better than the analytical model. The differences between the analytical and experimental responses can have several sources. The first and most important reason is the level of assumptions made formulation of the analytical model. A second reason is that the initial imperfections in the 3D printed cylinder are not those assumed in the analytical model.

### 4.2 Use of Semi-Analytical Model for Design.

This section aims to show how the developed analytical model can be used for the design of NSD cylinders with a desired post-buckling response under axial compression. This can be accomplished by following some steps of the framework in Fig. 4 backward. Maps for design parameters need to be developed first. To do so, numerous panels with varying geometry should be analyzed for their response to compression. This enables the selection of panels that fulfill the design goal. Of course, the selected panels from design maps represent the equivalent panels that were described in Sec. 2. It must be noted that the dimensions of the thickened areas on the top and bottom of the buckling panels affect some design factors, such as its stiffness, and the axial shortening corresponding to the buckling event. These complexities make the design a trial-and-error process.

#### 4.2.1 Design Maps for Post-Buckling Response of Cylindrical Panels.

The development of design maps for panels within a range of thickness, radius, and material properties is presented in this subsection. It is acknowledged that the range of panel height and width considered is limited. Providing generalized design maps, a detailed design guideline for the NSD cylinder requires a separate study that the authors plan to conduct in the future.

The analytical model for cylindrical panels presented in Sec. 3 was used to develop the design maps. A parametric evaluation was conducted first considering the post-buckling response of 16 cylindrical panels with changes in width and length from 20 mm to 80 mm. The radius and thickness of all panels are 50 mm and 0.5 mm, respectively. The material is ABS with elastic modulus (E) of 1300 N/mm2 and Poisson's ratio (ν) equal to 0.2. The simulation was conducted by considering the contribution of 100 coupled modes (Λ = ψ = 4).

The axial force-deformation responses for the cases outlined above are plotted in Fig. 9. A visual examination of the traces shows that the displacement corresponding to the first buckling event is higher for narrower panels. Also, for most of the cases, longer panel have a larger drop in stress due to buckling. Panels with width (b) or length (l) equal to 20 mm show very smooth buckling without a significant drop in axial stress during the buckling event. Following the response of the panels with l = 20 mm shows that buckling in those panels occurs at higher strains than other ones. For panels with b = 20 mm and l > 20 mm, the difference in the buckling strain is not significant. As expected, the initial stiffness for all panels is the same, and narrower panels have a larger post-buckling stiffness.

Fig. 9
Fig. 9
Close modal

From a design perspective, Fig. 9 shows that it is not possible to design an NSD cylinder with a negative post-buckling stiffness if the dimensions of the buckling targeted panels are in the range of this parametric evaluation. To experience the slightest secondary stiffness, the suggestion is to divide the cylinder to fewer numbers of segments which are wider and longer. However, if the designer aims to see fastest significant buckling in the NSD cylinder, one of the panels should have an equivalent length of 20 mm and a width of 80 mm.

Design contours for effective design parameters (Fig. 10) were developed using the information from Fig. 9. Contours are provided for four design parameters: (Fig. 10(a)) buckling strain, (Fig. 10(b)) buckling stress, (Fig. 10(c)) stress drop, and (Fig. 10(d)) secondary stiffness. The contour plots were generated by fitting a surface, for each design parameter, to the 16 data points from the parametric evaluation. The contour plots can help select the dimension of buckling panels based on desired design parameters.

Fig. 10
Fig. 10
Close modal

#### 4.2.2 Design of Non-Uniform Stiffness Cylinders With Targeted Post-Buckling Response.

This section illustrates the design of an NSD cylinder with a targeted elastic post-buckling response under uniform axial compression using the design maps in Fig. 10. The cylinder is to have an inner radius of 50 mm, a height of 100 mm, and the thickness of the buckling areas equal to 0.5 mm. The material is linear elastic with a modulus of 1300 N/mm2 and Poisson's ratio of 0.2. The cylinder is required to experience three local buckling events when the global axial shortening reaches 0.24 mm, 0.32 mm, and 0.4 mm. Also, b, l1, l2, and b1 (see Fig. 1) in all three segments are 60 mm, 30 mm, 10 mm, and 30 mm, respectively.

The displacement corresponding to the buckling in segment A (see Fig. 2) is calculated as follows:
$δ0A=ϵ0l1A′+2σ0btbk1A=0.24mm$
(18)

Here, $δ0A$ is the end shortening corresponding to buckling in segment A. Also, $l1A′$ and $k1A$ are expressed in Eqs. (2) and (3). In Eq. (18), ε0 is a function of $l(l1′A)$ and b; $l1A′$ is a function of l1, tb, and t1A; σ0 is a function of $l(l1′A)$ and b; and k1A is a function of b, b1, l1, l2, tb, and t1A. Thus, $δ0A$ is a function of six independent variables, from which five are predefined here and only t1A remains unknown.

By following a trial-and-error process, one can design the cylindrical panel that buckles under 0.24 mm axial shortening. However, this process is time-consuming. Thus, a numerical minimization problem using a genetic algorithm was formulated using the software mathematica [34]. The objective was to minimize ($δ0A$ − 0.24)2. The same process was used to find t1B and t1C by numerically minimizing ($δ0B$ − 0.32)2 and ($δ0C$ − 0.40)2. The obtained solutions were t1A = 3.89 mm, t1B = 1.30 mm, and t1C = 0.9 mm.

Segments D do not affect the design objective, and they just should be stiff enough to separate the buckling events in segments A, B, and C. A thickness of 1.0 mm for t2 and a length of 60 mm for l3 separated the buckling segments well in previous experiments (Sec. 4.1), and thus, the same values were used here.

A finite element model for the NSD cylinder with the designed dimensions was developed in the program abaqus to examine the accuracy of the design process. The modeling details were as described in Sec. 4.1.

The axial response of the designed NSD cylinder from the finite element analysis is shown in Fig. 11(a), and contour maps of the radial deformations are shown in Fig. 11(b). The buckling events happen in the designed order, and the axial shortening corresponding to the buckling events are well matched between the design target and the analysis results. This example shows that it is possible to design NSD cylinders with a targeted post-buckling response using the presented analytical model. However, further studies are necessary to develop the complete design procedure for any desired post-buckling response in NSD cylinders.

Fig. 11
Fig. 11
Close modal

## 5 Conclusions

Cylindrical shells with non-uniform stiffness distribution (NSD cylinders) are designed to exhibit localized buckling in targeted areas that are surrounded by thicker (and elastically responding) regions. In this study, the buckling events in targeted areas were assumed independent from each other and a semi-analytical model was developed to predict the post-buckling response of NSD cylinders under uniform axial compressive force. The model has three key elements: (i) subdivision of the cylinder into parallel panels with linear and nonlinear response, (ii) a semi-analytical approximate model to predict the nonlinear response of the cylindrical panels, and (iii) integrating the response from the individual panels. The axial post-buckling response of several NSD cylinders was predicted by the proposed model and compared with the axial response from experimental tests on 3D printed units and FEA. The comparison shows that the model presented in this paper accurately predicts the axial shortening corresponding to each buckling event in the cylinder. The results show that in some cases, the semi-analytical model predicted the post-buckling response of NSD cylinders better than the FEA, without special considerations in the FE modeling such as mesh effects, shear locking, hourglass problems, the effect of element type, and the influence of neglecting the panel thickness. However, it slightly underestimates the stress drop in the buckling events, and it slightly overestimates the post-buckling stiffness of the cylinders.

Finally, the functionality of this semi-analytical model for design was shown by developing response contours and using them to design an NSD cylinder with a desired post-buckling response. The post-buckling response of the NSD cylinder designed using the semi-analytical model was compared with the response of the cylinder from FEA. This comparison showed that the proposed semi-analytical model has the potential to design NSD cylinders for desired post-buckling responses with acceptable accuracy.

## Acknowledgment

The presented work was carried out with the support from the U.S. National Science Foundation under grant CMMI-1463164. The reported experiments were conducted by Mr. Jun Guo as part of a related research effort.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## References

1.
Fan
,
H.
,
2019
, “
Critical Buckling Load Prediction of Axially Compressed Cylindrical Shell Based on Non-Destructive Probing Method
,”
Thin-Walled Struct.
,
139
, pp.
91
104
.
2.
Emam
,
S. A.
, and
Inman
,
D. J.
,
2015
, “
A Review on Bistable Composite Laminates for Morphing and Energy Harvesting
,”
ASME Appl. Mech. Rev.
,
67
(
6
), p.
60803
.
3.
Wang
,
Y.
, and
Minor
,
M. A.
,
2018
, “
Design and Evaluation of a Soft Robotic Smart Shoe for Haptic Terrain Rendering
,”
IEEE/ASME Trans. Mechatron.
,
23
(
6
), pp.
2974
2979
.
4.
Roh
,
J.-H.
,
Oh
,
I.-K.
,
Yang
,
S.-M.
,
Han
,
J.-H.
, and
Lee
,
I.
,
2004
, “
Thermal Post-Buckling Analysis of Shape Memory Alloy Hybrid Composite Shell Panels
,”
Smart Mater. Struct.
,
13
(
6
), pp.
1337
1344
.
5.
Chang-qing
,
C.
, and
Ya-peng
,
S.
,
1997
, “
Stability Analysis of Piezoelectric Circular Cylindrical Shells
,”
ASME J. Appl. Mech.
,
64
(
4
), pp.
847
852
.
6.
Lee
,
H. J.
, and
Lee
,
J. J.
,
2000
, “
A Numerical Analysis of the Buckling and Postbuckling Behavior of Laminated Composite Shells With Embedded Shape Memory Alloy Wire Actuators
,”
Smart Mater. Struct.
,
9
(
6
), pp.
780
787
.
7.
Kundu
,
C. K.
,
Maiti
,
D. K.
, and
Sinha
,
P. K.
,
2007
, “
Post Buckling Analysis of Smart Laminated Doubly Curved Shells
,”
Compos. Struct.
,
81
(
3
), pp.
314
322
.
8.
Liu
,
S.
,
,
A. I.
, and
Burgueño
,
R.
,
2019
, “
Architected Materials for Tailorable Shear Behavior With Energy Dissipation
,”
Extrem. Mech. Lett.
,
28
, pp.
1
7
.
9.
Liu
,
D.
,
Kitipornchai
,
S.
,
Chen
,
W.
, and
Yang
,
J.
,
2018
, “
Three-Dimensional Buckling and Free Vibration Analyses of Initially Stressed Functionally Graded Graphene Reinforced Composite Cylindrical Shell
,”
Compos. Struct.
,
189
, pp.
560
569
.
10.
Liu
,
S.
,
,
A. I.
, and
Burgueño
,
R.
,
2018
, “
Energy Harvesting From Quasi-Static Deformations via Bilaterally Constrained Strips
,”
J. Intell. Mater. Syst. Struct.
,
29
(
18
), pp.
3572
3581
.
11.
Kobayashi
,
T.
,
Mihara
,
Y.
, and
Fujii
,
F.
,
2012
, “
Path-Tracing Analysis for Post-Buckling Process of Elastic Cylindrical Shells Under Axial Compression
,”
Thin-Walled Struct.
,
61
, pp.
180
187
.
12.
Wagner
,
H. N. R.
,
Hühne
,
C.
, and
Niemann
,
S.
,
2020
, “
Buckling of Launch-Vehicle Cylinders Under Axial Compression: A Comparison of Experimental and Numerical Knockdown Factors
,”
Thin-Walled Struct.
,
155
, p.
106931
.
13.
Wagner
,
H. N. R.
,
Hühne
,
C.
,
Niemann
,
S.
, and
Khakimova
,
R.
,
2017
, “
Robust Design Criterion for Axially Loaded Cylindrical Shells—Simulation and Validation
,”
Thin-Walled Struct.
,
115
, pp.
154
162
.
14.
Wagner
,
R.
,
2019
, “
Robust Design of Buckling Critical Thin-Walled Shell Structures
,”
Doctoral Dissertation
,
DLR Dtsch. Zent. fur Luft- und Raumfahrt e.V.—Forschungsberichte
.
15.
Degenhardt
,
R.
,
Castro
,
S. G. P.
,
Błachut
,
J.
,
Arbelo
,
M. A.
, and
Khakimova
,
R.
,
2017
, “Stability of Composite Shell-Type Structures,”
Stability and Vibrations of Thin-Walled Composite Structures
, pp.
253
428
.
16.
Castro
,
S. G. P.
,
Zimmermann
,
R.
,
Arbelo
,
M. A.
,
Khakimova
,
R.
,
Hilburger
,
M. W.
, and
Degenhardt
,
R.
,
2014
, “
Geometric Imperfections and Lower-Bound Methods Used to Calculate Knock-Down Factors for Axially Compressed Composite Cylindrical Shells
,”
Thin-Walled Struct.
,
74
, pp.
118
132
.
17.
Ning
,
X.
, and
Pellegrino
,
S.
,
2015
, “
Imperfection-Insensitive Axially Loaded Thin Cylindrical Shells
,”
Int. J. Solids Struct.
,
62
, pp.
39
51
.
18.
Cox
,
B. S.
,
Groh
,
R.
, and
Pirrera
,
A.
,
2019
, “
Nudging Axially Compressed Cylindrical Panels Towards Imperfection Insensitivity
,”
ASME J. Appl. Mech.
,
86
(
7
), p.
071010
.
19.
Hu
,
N.
, and
Burgueño
,
R.
,
2015
, “
Tailoring the Elastic Postbuckling Response of Cylindrical Shells: A Route for Exploiting Instabilities in Materials and Mechanical Systems
,”
Extrem. Mech. Lett.
,
4
, pp.
103
110
.
20.
Burgueño
,
R.
,
Hu
,
N.
,
Heeringa
,
A.
, and
Lajnef
,
N.
,
2014
, “
Tailoring the Elastic Postbuckling Response of Thin-Walled Cylindrical Composite Shells Under Axial Compression
,”
Thin-Walled Struct.
,
84
, pp.
14
25
.
21.
Hu
,
N.
, and
Burgueño
,
R.
,
2015
, “
Elastic Postbuckling Response of Axially-Loaded Cylindrical Shells With Seeded Geometric Imperfection Design
,”
Thin-Walled Struct.
,
96
, pp.
256
268
.
22.
Hu
,
N.
, and
Burgueño
,
R.
,
2017
, “
Harnessing Seeded Geometric Imperfection to Design Cylindrical Shells With Tunable Elastic Postbuckling Behavior
,”
ASME J. Appl. Mech.
,
84
(
1
), p.
11003
.
23.
Hu
,
N.
, and
Burgueño
,
R.
,
2018
, “
Cylindrical Shells With Tunable Postbuckling Features Through Non-Uniform Patterned Thickening Patches
,”
Int. J. Struct. Stab. Dyn.
,
18
(
02
), p.
1850026
.
24.
Guo
,
J.
,
Liu
,
S.
, and
Burgueño
,
R.
,
2017
, “
Tailoring the Elastic Postbuckling Response of Thin-Walled Cylindrical Shells for Applications in Mechanical Devices and Adaptive Structures
,”
Proceedings of the ASME 2017 Conference on Smart Materials, Adaptive Structures and Intelligent Systems. Volume 2: Modeling, Simulation and Control of Adaptive Systems; Integrated System Design and Implementation; Structural Health Monitoring
,
Snowbird, UT
,
Sept. 18–20
, p.
V002T03A037
.
25.
Koiter
,
W. T.
,
Elishakoff
,
I.
,
Li
,
Y. W.
, and
Starnes
,
J. H.
, Jr.
,
1994
, “
Buckling of an Axially Compressed Cylindrical Shell of Variable Thickness
,”
Int. J. Solids Struct.
,
31
(
6
), pp.
797
805
.
26.
Chen
,
L.
,
Rotter
,
J. M.
, and
Doerich
,
C.
,
2011
, “
Buckling of Cylindrical Shells With Stepwise Variable Wall Thickness Under Uniform External Pressure
,”
Eng. Struct.
,
33
(
12
), pp.
3570
3578
.
27.
Chen
,
Z.
,
Yang
,
L.
,
Cao
,
G.
, and
Guo
,
W.
,
2012
, “
Buckling of the Axially Compressed Cylindrical Shells With Arbitrary Axisymmetric Thickness Variation
,”
Thin-Walled Struct.
,
60
, pp.
38
45
.
28.
Brar
,
G. S.
,
Hari
,
Y.
, and
Williams
,
D. K.
,
2008
, “
Buckling of Axisymmetric Cylindrical Shells of Variable Thickness: Finite Difference Solution
,”
Proceedings of the ASME 2007 Pressure Vessels and Piping Conference. Volume 1: Codes and Standards
,
San Antonio, TX
,
July 22–26
,
American Society of Mechanical Engineers
, pp.
627
633
.
29.
Fan
,
H. G.
,
Chen
,
Z. P.
,
Feng
,
W. Z.
,
Zhou
,
F.
,
Shen
,
X. L.
, and
Cao
,
G. W.
,
2015
, “
Buckling of Axial Compressed Cylindrical Shells With Stepwise Variable Thickness
,”
Struct. Eng. Mech.
,
54
(
1
), pp.
87
103
.
30.
Greiner
,
R.
,
1972
,
Ein baustatisches Lösungsverfahren zur Beulberechnung dünnwandiger Kreiszylinderschalen unter Manteldruck, Bauingenieur-Praxis 17
,
Berlin
:
Ernst & Sohn Verlag
.
31.
Yamaki
,
N.
,
1984
,
Elastic Stability of Circular Cylindrical Shells
, Vol.
27
,
Elsevier
,
New York
.
32.
Shen
,
H. S.
,
2007
, “
Postbuckling Analysis of Axially Loaded Piezolaminated Cylindrical Panels With Temperature Dependent Properties
,”
Compos. Struct.
,
79
(
3
), pp.
390
403
.
33.
Simulia
,
2019
, “
Abaqus User’s Manual Version 2019
,”
Dassault Systèmes Simulia Corp.
,
Providence, RI
.
34.
I. Wolfram Research
,
2020
,
Mathematica
,
Version 12
,
Champaign, IL
.