## 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 [2–10]. 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 [12–15]. 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 [19–23]. 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.

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 [27–29]. 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 *t _{b}* in Fig. 1), where

*thicker*zones (areas with thickness equal to

*t*

_{2}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.

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 (*t _{b}* in Fig. 3) that yields the same deformation Δ under equal bending moment

*M*.

*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\u2032$ is the equivalent length of the same area, and

*t*and $t1i$ are the thicknesses of the panel's mid-section and the thicker areas,

_{b}*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):

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.

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

*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

*u*and

*v*are related to out-of-plane deformation

*w*and

*w**and the airy stress function if the following way:

*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):

Here, *x* and *y* are the in-plane coordinates of the panel, *i*_{1}, *i*_{2}, *j*_{1}, and *j*_{2} are the number of the waves in the *x*- and *y*-directions (respectively), *a _{i}*

_{1i2j1j2}is the unknown coefficient for the radial deformation,

*μ*

_{i}_{1i2j1j2}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.

*w*and

*w**from Eq. (9), Eq. (5) can be written as follows:

*ψ*are the maximum number of contributed modes in Eq. (9) for the

*x*- and

*y*-directions, and

*f*is the multiplier for the cosine function in ∇

_{ijk}^{4}

*f*. Solving Eq. (10) results in a function

*F*that has similar cosine functions to ∇

^{4}

*f*but different multipliers to those cosine expressions (

*F*)

_{ijk}*λ*(

*y*) and

*θ*(

*x*) with ∇

^{4}

*λ*(

*y*) = ∇

^{4}

*θ*(

*x*) = 0. Then, it can be expressed as follows:

*N*over width of the panel should be equal to the applied force (

_{x}*P*). Thus,

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

*θ*(

*x*) can be expressed as

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 *a _{i}*

_{1i2j1j2}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/mm^{2} 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.

Cases | Length (l) | Width (b) | Radius (R) | Thickness (t) |
---|---|---|---|---|

mm | mm | mm | mm | |

(a) | 80 | 80 | 50 | 0.5 |

(b) | 80 | 40 | 50 | 0.5 |

(c) | 40 | 60 | 50 | 0.5 |

(d) | 60 | 40 | 50 | 0.5 |

Cases | Length (l) | Width (b) | Radius (R) | Thickness (t) |
---|---|---|---|---|

mm | mm | mm | mm | |

(a) | 80 | 80 | 50 | 0.5 |

(b) | 80 | 40 | 50 | 0.5 |

(c) | 40 | 60 | 50 | 0.5 |

(d) | 60 | 40 | 50 | 0.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 *μ _{i}*

_{1i2j1j2}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.015

*h*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.

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.

Test unit | t (mm)_{b} | t_{1A} (mm) | t_{1B} (mm) | t_{1C} (mm) | t_{2} (mm) | l_{1} (mm) | $l1A\u2032$ (mm) | $l1B\u2032$ (mm) | $l1C\u2032$ (mm) | l_{3} (mm) | b_{1} (mm) | b_{2} (mm) |
---|---|---|---|---|---|---|---|---|---|---|---|---|

CYL-1 | 0.5 | 1 | 1.3 | 4 | 1 | 30 | 10.56 | 7.15 | 1.33 | 60 | 35 | 65 |

CYL-2 | 0.5 | 1 | 1.5 | 4 | 1 | 30 | 10.56 | 5.75 | 1.33 | 60 | 35 | 65 |

CYL-3 | 0.5 | 1 | 2 | 4 | 1 | 30 | 10.56 | 3.75 | 1.33 | 60 | 35 | 65 |

CYL-4 | 0.5 | 1 | 2.5 | 4 | 1 | 30 | 10.56 | 2.65 | 1.33 | 60 | 35 | 65 |

Test unit | t (mm)_{b} | t_{1A} (mm) | t_{1B} (mm) | t_{1C} (mm) | t_{2} (mm) | l_{1} (mm) | $l1A\u2032$ (mm) | $l1B\u2032$ (mm) | $l1C\u2032$ (mm) | l_{3} (mm) | b_{1} (mm) | b_{2} (mm) |
---|---|---|---|---|---|---|---|---|---|---|---|---|

CYL-1 | 0.5 | 1 | 1.3 | 4 | 1 | 30 | 10.56 | 7.15 | 1.33 | 60 | 35 | 65 |

CYL-2 | 0.5 | 1 | 1.5 | 4 | 1 | 30 | 10.56 | 5.75 | 1.33 | 60 | 35 | 65 |

CYL-3 | 0.5 | 1 | 2 | 4 | 1 | 30 | 10.56 | 3.75 | 1.33 | 60 | 35 | 65 |

CYL-4 | 0.5 | 1 | 2.5 | 4 | 1 | 30 | 10.56 | 2.65 | 1.33 | 60 | 35 | 65 |

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.

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.05*t*. 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/mm^{2}.

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.

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/mm^{2} 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.

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.

#### 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/mm^{2} 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*, *l*_{1}, *l*_{2}, and *b*_{1} (see Fig. 1) in all three segments are 60 mm, 30 mm, 10 mm, and 30 mm, respectively.

Here, $\delta 0A$ is the end shortening corresponding to buckling in segment A. Also, $l1A\u2032$ and $k1A$ are expressed in Eqs. (2) and (3). In Eq. (18), *ε*_{0} is a function of $l(l1\u2032A)$ and *b*; $l1A\u2032$ is a function of *l*_{1}, *t _{b}*, and

*t*

_{1A};

*σ*

_{0}is a function of $l(l1\u2032A)$ and

*b*; and

*k*

_{1A}is a function of

*b*,

*b*

_{1},

*l*

_{1},

*l*

_{2},

*t*, and

_{b}*t*

_{1A}. Thus, $\delta 0A$ is a function of six independent variables, from which five are predefined here and only

*t*

_{1A}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 ($\delta 0A$ − 0.24)^{2}. The same process was used to find *t*_{1B} and *t*_{1C} by numerically minimizing ($\delta 0B$ − 0.32)^{2} and ($\delta 0C$ − 0.40)^{2}. The obtained solutions were *t*_{1A} = 3.89 mm, *t*_{1B} = 1.30 mm, and *t*_{1C} = 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 *t*_{2} and a length of 60 mm for *l*_{3} 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.

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