Bending Stability of Corrugated Tubes With Anisotropic Frustum Shells

Snap-through instabilities of elastic shells offer a novel method to create reconfigurable structures and geometrically adaptable systems. Recent studies have explored the underlying mechanics of different elastic shell geometries [1–4], and research has also harnessed their instabilities for different morphological applications [5,6]. Among common geometries, the open-top conical frustum shell can be inverted to generate non-trivial stable configurations [7]. In their most popular implementation, open-top frusta are arranged in an opposing fashion to create a functional corrugated tube commonly known as the bendy or flexible drinking straw (Fig. 1). This corrugated system allows for morph-and-lock features with both axial and rotational deformations. A full inversion of the frustum can axially shorten the corrugated tube while a partial inversion can allow for global rotation. The separate frusta can undergo different inversions, allowing for versatile shape morphing over the length of the corrugated tube.

Corrugated tubes made of open-top frusta have found practical applications beyond the simple drinking straw. The same principles are used to create bain circuits [8,9], hoses for transporting liquids and gases [10,11], and covers for moving equipment [12]. The ability to elongate, shorten, and adaptively bend allows this subset of corrugated tubes to conform to different shapes, move into new configurations, and to connect points in otherwise hard to navigate environments (e.g., within the human body or a car engine). Furthermore, the corrugations increase the stiffness of the tubes for external out-of-plane loads and increase the capacity for the tubes to resist orthogonal flattening or crushing. Recently, these corrugated tubes of open-top frusta have also been investigated for novel applications in engineering. The axial bi-stability of the multiple successive frusta in the tube has been used to create energy absorbing and adaptable metamaterials [13,14], and the bending adaptivity has been proposed for use in creating versatile origami-inspired finger grippers and crawling robots [15]. Bending bi-stability of the frusta allows the tube to be reconfigured to arbitrary 3D space curves, which opens the door for transdimensional deformations [16].

Despite the long history of practical applications with bending straws [17], the underlying mechanics of these systems have not been explored in much detail, especially with respect to the bending stability. The work by Zhang explored partial inversions of single open- and closed-top frusta under lateral point loads [7]. The work evaluated critical forces and the bi-stable transitions through experiments and a series of parametric studies on the geometry of the frusta. A different study by Bende et al. explored the bending stability of two opposing frusta (a unit cell of a corrugated tube) under a downward load placed non-axially at the edge of the frusta [18]. For a limited set of geometries, the work showed that internal pre-stress or overcurvature is necessary to achieve a stable bent state. The previous studies, however, have had several limitations. They have only explored a limited set of geometries and have not evaluated the influence of the “crease” which connects two separate conical frusta. From the experimental perspective, an extensive parametric study has not been feasible for the large number of geometric parameters, as it would require a large number of samples and testing time. In the previous studies, the partial inversions were activated by non-axial point loads to investigate the stability behavior of a single frustum rather than of a corrugated tube consisting of multiple frusta. While this loading mode was easy to achieve experimentally, it differs from the practical use of flexible bending straws which would typically experience bending moments. Finally, while the properties of circular frusta have been partially explored, cross-sectional influences on the bending stability have not been studied in previous work. Our work in this paper will give insight to these additional topics of interest.

The paper is organized as follows. Section 2 describes the geometry of the corrugated tubes, the adapted bar-and-hinge method used for simulations, and the four-point bending set-up used to explore bending bi-stability. In Sec. 3, the bending stabilities of circular cross sections are explored, and a simple mechanism is used to explain the underlying mechanics of the system. Section 4 evaluates the bending stability of modified cross sections and gives insight to how cross-sectional anisotropy, concavity, and non-symmetry can influence the stable configurations and energy landscape. Physical models with orthotropic bending stabilities are presented in Sec. 5. A discussion and conclusion addressing the findings, limitations, and potential extensions is presented in Sec. 6.

A typical instance of the bendy straw is composed of repeating units, where each unit cell consists of two different frusta that are connected by a crease line (Fig. 1(a)). In certain instances, the unit cell can snap between two axial stable states (Fig. 1(b)), extended and retracted, through a full inversion of the top frustum [13,18]. This type of axial bi-stability allows a straw to deploy and to maintain configurations at multiple lengths [13]. Another type of possible bi-stable behavior is through a partial inversion of the top frustum, which can provide continuously varied bending configurations for the straw, and here is referred to as bending stability [18].

The shape of the straw will depend on the longitudinal geometry of the unit cell (Fig. 1(c)) and the shape of the tube cross section. The longitudinal shape can be defined by the following five parameters: slant angles of (1) the top (θ₁) and (2) the bottom (θ₂) frusta, the radius of (3) the outer (R) and (4) the inner (r) edge, and (5) the shell thickness (t). The slope of the bottom frustum θ₂ is greater than θ₁ to avoid self-intersection during full inversion of the top frustum in axial compression [18]. Futhermore, r must be identical for all units in a multi-cell straw (Fig. 1(a)) because each unit must be connected to its adjacent unit via the inner edge. While in principle the other four parameters can be varied independently for each unit cell, we assume they are constant for the entire corrugated tube. To constrain our study to classical straw structures similar to those studied in Refs. [7,13,18], we set the ratio of the outer to inner radius to be R/r = 1.15, and the angle of the bottom frustum to be θ₂ = 70 deg. The outer radius is set to 31 mm. We then normalize the thickness with the outer radius as t/R. The normalized thickness (t/R) and the slant angle of the top frustum are systematically varied to explore the geometric influence on stability.

The parameter λ_φ serves a similar role to “a radius of curvature” in defining the negatively curved region in the dumbbell shape. The shape becomes a track shape when λ_φ → ∞, and the curvature of the negatively curved region increases as λ_φ decreases. The equation for the outer edge is formulated by replacing r with R in Eq. (2), and a consistent gap of R − r is created. The definitions in Eqs. (2) and (3) are chosen such that we obtain a smooth transition from the negatively curved region to the adjacent semi-circles. This shape is named as dumbbell shape because of its similarity to a physical dumbbell.

The bar and hinge model is a reduced-order simulation tool that has been used to approximate the loading response of thin-sheet structures, based on their global geometry and material properties [19,20]. The model can provide a rapid approximation of nonlinear mechanical behaviors with much less convergence issues than with conventional finite element models (discussed later in this section). The model was originally developed for analysis of straight-crease origami structures [21] and has been recently modified to simulate arbitrary sheets and origami with curved creases [22,23]. The bar and hinge model was carefully examined under various loading modes, and the predictions are in good agreement with those of theory, FE simulation, and experiments [22,23]. Curved crease origami is geometrically and physically similar to the open top frusta explored in this work. Therefore, the model is employed as the main tool to investigate the bending stability. In the bar and hinge model, the thin-sheet structure is discretized into nodes, bars, and hinges (Fig. 2(a)). There are three different elements for the main deformation behaviors: (1) straight bars that capture in-plane stretching and shearing, (2) bending hinges that capture out-of-plane bending of the sheet, and (3) folding hinges that capture rotations of the creases that connect separate sheets. The stiffness of each of the above three elements is derived based on linear-elastic material properties and mesh geometry [23]. Here, we use Y = 1300 MPa, ν = 0.45 for Young’s modulus and Poisson’s ratio, based on typical properties of isotropic and homogeneous polypropylene (PP) [7]. The straw geometry is discretized into triangular segments based on a metric called the aspect ratio, α, which describes the length to width ratio of the segments and thus controls the density of the mesh [23]. For each triangular segment, the aspect ratio α is defined as the ratio of the side length perpendicular to the fold, H, to the side length roughly parallel to the fold, W (Fig. 2(a)). A larger aspect ratio α results in a finer mesh discretization. As compared to the finite element method, it takes much less time to run simulations with the bar and hinge model, due to the fewer degrees-of-freedom. Prior work using the bar and hinge model to simulate thin sheet structures has successfully approximated the stiffness and deformation behaviors of different corrugated sheet structures [22,23].

Complicated behaviors affecting the folding stiffness are simplified into a single, case-specific parameter, called the length scale factor L*. This factor is found to be a linear function of the sheet thickness [25]. Given that the thickness t is normalized by R for a scalable analysis without consideration of the folds, L* is also normalized as L*/R for a fully scalable analysis with respect to folding stiffness. Unless mentioned specifically, in this work we define the folding hinges to have a stiffness of zero (L*/R = ∞) to focus the research on the geometrical influence on the stability behaviors of the thin-shell frusta.

We use a four-point bending test to investigate the bending response of the corrugated tube, and a two-unit straw is chosen because it is the simplest system that allows interactions within units. The potential boundary effect with respect to the number of the units will be discussed in the next subsection. Two stiffer hollow sections are connected to the ends of the straw via creases (Fig. 2(c), dark area). These sections provide a span between the loading points and the supports. These end sections are set to be 10⁶ stiffer than the straw so their deformations are negligible in comparison. The two nodes on the left end of the tube that lie on the x − z plane are restrained in all translational directions (Fig. 2(c)), while the two nodes on the right end (also on the x − z plane) are only restrained in x and y to allow global longitudinal movement of the tube. A displacement controlled simulation is performed by prescribing a downward y − displacement to the four nodes that connect the frusta to the hollow sections and also intersect the x − z plane (Fig. 2(c)). Bending about a different axis is realized by rotating the structure about the z-axis, without changing the criterion of assigning boundary and loading conditions. Total energy is calculated for every increment of the displacement controlled simulation by summing the strain energy in the bars and the hinges. The rotational angle ϕ is computed for every increment.

With the bar and hinge model, the structural stiffness can be overestimated for cases where torsional deformations are present [23]. Because torsion is involved in the bending of the straw corrugation (Fig. 2(c)), a calibration of the model is necessary to obtain reliable results. A finite element (FE) model (abaqus/standard [26]) of the two-unit circular straw is built using shell elements and rotational hinges to provide high-fidelity simulation results for the calibration (Fig. 2(d)). Shells are meshed with S4 general purpose elements, where adjacent frusta are connected via connector elements with prescribed rotational stiffness to simulate crease rotations (the same stiffness as defined in Eq. (5)). Each end of the corrugated portion is connected to a cylindrical tube with 10⁶ higher Young’s modulus. In the FE simulation, the bending is simulated by applying the same four-point bending setup and the loading scheme is controlled by the modified Riks method [27]. Linear-elastic material properties are set to be the same as the bar and hinge model. The rotational angle ϕ and the strain energy are computed in every iterations. Convergence with respect to the mesh density is examined, and a mesh size of roughly 0.03R × 0.03R is chosen because it provides a strain energy prediction that is within 0.2% of a mesh with 0.015R × 0.015R elements.

Because torsional deformations are expected in the simulations (Fig. 2(c)), an aspect ratio α ≈ 1 is used in the mesh of the bar and hinge model as is recommended in Ref. [23]. To calibrate the bar and hinge model, the energy profile (Fig. 2(e)) is compared to the FE results to quantify the difference. In order to be consistent with the scalable analysis, the strain energy E is non-dimensionalized as $E¯=106⋅E/(Y⋅V3)$⁠, where Y is the Young’s modulus and V denotes the material cost of the corrugated part. The factor of 10⁶ is included to make the figures that display non-dimensionalized energy more clear. For a bi-stable bending process, the energy profile has an energy barrier $E¯b$ and the valley point that corresponds to the stable rotation ϕ_b. Setting the discrepancy of energy barriers $E¯b$ and the secondary stable state ϕ_b as the minimizing objectives, the calibration returns stretching and bending stiffness as shown in Eqs. (7) and (8).

With the calibration, a close match of energy barriers $E¯b$ (Fig. 2(f)) and the second stable state ϕ_b (Fig. 2(g)) are achieved for θ₁’s in the range from 15 deg to 27.5 deg and thicknesses t/R = 0.013, 0.0145, 0.016, 0.0175. The largest relative error of the bar and hinge model for predicting energy barriers when compared to the FE results is 3.8%. The second stable state ϕ_b shows a nearly linear relationship versus θ₁, where the largest error is reported as 3.3%. The energy barrier $E¯b$ increases with the θ₁ (Fig. 2(f)), since more energy input is required for partially inverting a larger frustum. A corrugation with larger θ₁ achieves the second stable state at a larger magnitude of rotation. The top frustum gets taller and a larger rotation is required to trigger buckling and inversion.

In order to examine the bar and hinge model’s accuracy for predicting bending stability of systems with different geometries, 25 cases with various θ₁ and t/R are simulated and compared with the FE results. The stability results are presented in Fig. 2(h), where the color in each grid shows the stability type of the structure (i.e. mono-stable or bi-stable). Bending stabilities of 92% of the different geometric cases are predicted to be the same with the bar and hinge model and the FE models. Although two cases near the transition boundary between mono- and bi-stability are characterized differently with the two models, the global trends for the stability behavior is the same. This stability map will be discussed in the next section, regarding the stability trends with respect to θ₁ and shell thickness. Simulating the 25 cases with FE took ∼18 h, while the same simulations with the bar and hinge model took only 84 minutes. This improved computational efficiency with the bar and hinge model is necessary for the parametric variations in the following sections where we perform hundreds of variations varying the geometries of the systems. The model is still capable of accurately estimating the global behavior of the systems.

In the four-point bending simulation, each end of the corrugated tube is connected to a much stiffer hollow cylinder to provide a span between the loading points and the supports. Although the deformations and strain energies of these hollow sections are negligible as shown in later figures, they introduce stiff constraints on the conical frusta that may alter the loading response. Therefore, corrugated tubes with more unit cells are tested (Fig. 3(a)), and a convergence study on number of units is employed to explore the boundary effect. Here, a force-controlled simulation is implemented using the Modified Generalized Displacement Control Method (MGDCM) [28]. This algorithm is constructed based on the arc length method, and it can accurately follow the equilibrium path after the tangent stiffness becomes negative.

For the four-point bending test with a force control, there will be a nearly flat moment stage between the loading points, with the center moment value denoted by M (Fig. 3(b)). When the tube deforms, the frusta will snap and the moment value will fluctuate. The relationship between M and the tube rotation ϕ is illustrated in (Fig. 3(c)) where mono-stable systems will only have positive moments. The moment value of first limit point M_p is selected to compare the behavior of tubes with different numbers of units. Comparing the energy would be misleading as tubes with more units will experience more deformations and will contain more energy.

To study the convergence behavior, M_p is found for systems with different number of units, and a relative error of M_p is computed with respect to a system with N = 10 units. The results show a downward trend in error with respect to the number of units N (Fig. 3(d), note the logarithmic y-scale). When there is only one unit in the tube (N = 1), one of the nearly rigid cylinders is directly connected to the inverting frustum, and it undergoes non-negligible deformation when the frustum snaps. The substantial deformation due to the constraint will magnify the bending moment, and the relative error is near 100% for four combinations of different θ₁ and t/R. When there are two or more units in the tube, the relative error of M_p under various geometries is less than 3% (Fig. 3(d), see inset). Therefore, the boundary effect from the stiff constraints will substantially change the bending response when N = 1, but it can be reasonably neglected for N ⩾ 2. In other sections of this work, the stability of the system is explored using corrugated tubes with N = 2, based on the balance of accuracy and computational complexity.

Finally, the stability predictions for the same twenty five cases as in Fig. 2(h) are shown for tube with N = 10 in Fig. 3(e). For all but two cases on the boundary, the stability type is consistently predicted, regardless of the number of units.

This section focuses on how the bending stability of a corrugation with a circular cross section can be affected by geometric parameters defining the tube and the stiffness of the crease connecting separate frusta. Figure 2(c) shows that partial inversion of the frustum of a unit cell can be triggered with a certain magnitude of prescribed displacement. The top frustum is heavily deformed, whereas the shape of the bottom frustum stays relatively unchanged, suggest that bending responses can be affected by the geometry of the top frustum (see Figs. 1(b) and 1(c)), i.e., the slant angle θ₁. Given a two-unit structure with a circular cross section (θ₁ = 22.5 deg, t/R = 0.0175), the second stable state is achieved when the tube is rotated by ϕ_b ≈ 3.6 deg (Fig. 4(a), part 1). The total energy decreases after it reaches the energy barrier $E¯b$⁠, indicating the existence of a second stable state. The bending energy (local bending of the frusta) increases monotonically through the process, while the stretching energy shows a “snap-through” profile, suggesting that some stretching energy is released during the partial inversion. Thus, the stretching of the frusta is essential in having a bi-stable system. The secondary stable state disappears for a smaller θ₁ (Fig. 4(a), part 2). In this case, the release of stretching energy becomes smaller than the increase in bending energy, and there is no longer a snap-through.

The effects of θ₁ can be better understood by exploring how energy is redistributed during the partial inversion process. For case 1 and case 2 in Fig. 4(a), the change of energy between the peak point to the valley point of the stretching energy curve are shown in Fig. 4(b). The side shown in Fig. 4(b) is the compression side of the corrugated straw. A straw with a larger θ₁ has a longer and steeper conical frustum when other parameters are fixed. Compared to the frustum with smaller θ₁, such a frustum is heavily deformed, and the panels are stretched to a greater extent. Therefore, the case with θ₁ = 22.5 deg releases much more stretching energy (−2.53 mJ) around the partially-inverted frustum, compared to the case with θ₁ = 15 deg (only −0.44 mJ). In both cases, the redistributions show global increase of bending energy. While the larger θ₁ induces more gain of bending energy (Fig. 4(b)), the difference of the bending energy gain (2.27 mJ versus 0.86 mJ) is not big enough to counteract the release of stretching energy, so the case with larger θ₁ is bi-stable.

Another parameter that can affect bending stability is the shell thickness. The ratio of stretching to bending energy decreases as we increase the sheet thickness, because the bending energy scales with t³, while the stretching energy is proportional to t. Furthermore, stretching and bending have counteracting effects on the stability, where only stretching energy contributes to the bi-stable behavior. Thus, when we increase the thickness from t/R = 0.0175 to 0.0205, the total energy becomes monotonically increasing (Fig. 4(a), part 3 versus part 1). The stability of straws with other combinations of angle and thickness are shown in (Fig. 4(c)) and are characterized as bi-stable and mono-stable. Twenty-five cases of straws were investigated, where the frustum angle ranges as θ₁ ∈ {15 deg, 17.5 deg, 20 deg, 22.5 deg, 25 deg}, and the non-dimensionalized thickness ranges as t/R ∈ {0.0145, 0.016, 0.0175, 0.019, 0.0205}. The sheet thickness is normalized by the outer radius R, to make the analysis scalable. The bi-stable domain grows with an increase of θ₁ and shrinks with an increase of t, implying that a thinner corrugation with larger slant angle θ₁ favors bi-stable bending.

The stiffness of the crease lines which connect separate frusta can also affect the bending stability of the system. In part 4 of Fig. 4(a), a normalized crease stiffness of L*/R = 30 is added to the structure with θ₁ = 22.5 deg, and t/R = 0.0175 (equivalent geometry to structure with no folding stiffness in Fig. 4(a), part 1). The energy associated with folding of the crease increases monotonically during the partial inversion process and makes the system mono-stable. The effect of different crease stiffness on the bending stability of structures with different geometries is shown in Fig. 4(d). These expanded stability maps show one hundred combinations of θ₁ and t/R with three different length scale factors that define the folding stiffness in Eq. (5): L*/R = ∞, 30, 3. The folding stiffness is zero when L*/R = ∞, and the corresponding bi-stable domain reaches the maximum size. The bi-stable domain shrinks as L*/R decreases to 30 (meaning a finite folding stiffness). Moreover, only four cases are bi-stable when L*/R = 3 (a larger folding stiffness). Thus, for a corrugation with a reasonable folding stiffness, bending can only be bi-stable with significantly thinner shells and larger slant angles θ₁. In other words, crease stiffness in the corrugation resists bi-stable bending. These results also support the findings on limited geometries in Ref. [18]. These systems had a finite crease stiffness and only obtained bending stability due to pre-stress. Here, we show that with sufficient low folding crease stiffness, bending bi-stability is still possible by tuning the geometry.

In this section, we will verify the bending stability behaviors with a simplified mechanism model. The bar and hinge and FE model found that the following features affect bending stability. Bending bi-stability is only possible with partial inversion of the top frustum. The top frustum occupies a larger portion of the entire structure when θ₁ increases, and there is a higher percentage of released energy. Stretching energy is the driving component for a bi-stable strain energy profile, and it is proportional to the shell thickness t. In contrast, bending energy increases monotonically during the partial inversion, detracting from possible bi-stability, and the bending energy scales with t³. Therefore, decreasing the thickness t results in a second stable configuration to be possible for bending. Finally, the contribution of the creases is similar to that of the sheet bending, and creases with higher stiffness will result in less bi-stable structures.

Figure 5(a) shows the simplified planar four-bar linkage model at its initial state. The model consists of a horizontal rigid bar of length 2r, and two crank links with initial angles of α₁ = θ₁ and α₂ = π − θ₁. This frustum bending is modeled with two torsional springs (stiffness k₁) at the left and right ground supports (Fig. 5(a)). The springs are at a strain-free angle at the initial state, and their stiffness will be formulated to scale with t³, to represent relative scaling of bending of the thin sheet. Ground supports are distanced from each other by 2R, corresponding to the span of the outer edge (Fig. 1(c)). Our model is inspired by that shown in Ref. [18], however in our implementation, two translational springs are introduced to capture the stretching energy during the inversion. These springs are grounded via roller supports, and their top ends are connected to the top bar and cranks together. Stiffness of the translational spring is denoted by k₂ that will be formulated as a linear function of the shell thickness. The translational springs are strain-free at the length given by the initial configuration of the linkage model. As the mechanism approaches the inversion, the lengths of the springs go to zero (maximum strain energy), and after inversion, the length of the springs increase again (releasing strain energy). The initial state of the full mechanism is a stress-free stable configuration, corresponding to the extended state of a unit cell (Fig. 1(b)). A bent state through partial inversion of the crank links (Fig. 5(a)) refers to the bent state of the straw structure in Fig. 1(b). However, in the bent state, the springs will contain some strain energy; thus, the stability needs to be confirmed.

The relationship of E versus ψ can thus be plotted to investigate the stability of the system, as well as characteristic values including energy barriers $E¯b$ and the stable state ψ_b (Fig. 5(b)). The parameter space of the linkage model is {R, r, k₁, k₂}, and θ₁ as the initial value of α₁. The ratio for R/r is the same as with the other analyses, and θ₁ is varied in a broader range, from 15 deg to 55 deg. The stiffness of the translational springs is set to be k₂ = kt as a linear function of the shell thickness, where k takes into account the material property and the frustum dimensions that are irrelevant to t. While the choice of k affects the energy magnitude, it does not affect the stability prediction, the bi-stable configuration, nor the trend of the energy barrier, thus k can be arbitrary number for our formulation. The torsional stiffness k₁ is determined by k₂ = k₁/(ξ · t²), where the t² compliments dimensional difference between bending energy and stretching energy. The coefficient ξ determines the relative stiffness of stretching to bending.

The energy is non-dimensionalized as $E¯=E/(kt3)$⁠, so the stiffness parameter k can be eliminated from the energy expression. A typical bi-stable energy profile with tuned stiffness coefficient is shown in Fig. 5(b). The energy barrier $E¯b$ and the stable configuration ψ_b are obtained accordingly. The coefficient ξ is chosen to be 75 to match the transition boundary in the stability map (Fig. 5(c)) to that of the bar and hinge model. While the stability map does not match Fig. 4(c) exactly, the trends of stability with respect to θ₁ and t are captured consistently. The second stable state is possible with a higher slant angle θ₁ or a lower shell thickness t.

On the side where the top bar moves down (the left side in Fig. 5(a)), the crank is inverted during the bending process. The translational spring is first compressed to zero length and then released to a longer length afterwards. Correspondingly, the energy in the translational spring first increases then decreases, undergoing a bi-stable transition (Fig. 5(d)). On the same side, the torsional spring is rotated as the slope ψ increases, resulting in a monotonic increase in the bending energy within the system. The other translational and torsional springs gain and release a minor amount of energy compared to those on the inversion side (Fig. 5(d)). Therefore, the total bending energy is monotonically increasing, and stretching is the driving factor behind the bi-stable profile. Furthermore, these results shown that the compression side of the frustum dominates the behavior and determines whether the structure will be bistable. The resulting behavior matches that observed from the FE and bar and hinge simulations, and the findings are further backed up in Sec. 4.1.

The influence of θ₁ can also be better investigated with the simplified model. The length of the translational springs is W · sinθ₁ = (R − r) · tan θ₁, which increases faster than the traveling range of the torsional springs (proportional to the θ₁). Therefore, with larger θ₁, there is a larger contribution from the translational springs, which increase the relative amount of stretching energy and help achieve the bent stable state. The effect that a lower thickness increases bending bi-stability can also be observed with the simple mechanism model. As the thickness t is reduced the stretching energy (which is bi-stable) has a larger relative contribution to the total energy than the bending energy (which increases monotonically). The effect of thickness is also evident from Eq. (10), since k₁/k₂ grows quadratically with the thickness (k₁/(k₂ · t²) is set to be a constant).

The energy barrier $E¯b$ and stable state ψ_b both increase with θ₁ (Figs. 5(e), 5(f)), matching the observation from the bar and hinge model (Figs. 2(f), 2(g)). The linkage becomes taller when θ₁ increases while the width of the top bar is fixed, requiring a larger rotation to reach the stable configuration (Fig. 5(e)). Furthermore, the taller linkage generates longer translational springs, pushing up the energy barriers $E¯b$ (Fig. 5(f)). However, these changes are negligible compared to the effect of θ₁. The stable state ψ_b is smaller than those reported for the bar and hinge model (Fig. 5(e)), because the linkage only represents a single inverted frustum from the two-units corrugation (Fig. 2(c)). The energy barrier cannot be directly compared with the bar and hinge model because the arbitrary parameter k includes geometric effects not captured by the linkage, and the non-diminsionalization also needs to be different for the two models. Nonetheless, the energy and stable equilibria trends for θ₁ and t are captured consistently by the linkage model.

In order to analyze the effect of the corrugation cross section on the bending stability, three different variations of cross-sectional shapes are proposed in Sec. 2.1. These include a track shape where the circular cross section is elongated, a further modification to include sides with negative curvature, and an asymmetric cross section resembling a three-leaf clover. The first cross section variation that we investigated is the track shape, created by pulling two semi-circles away from the centroid and re-connecting them with straight lines. The parameter d is the distance between each semi-circle and the centroid. Starting with d = 0 (a circular cross section), an orthotropic cross section can be formed by adopting d > 0. The structure can then be bent about the two perpendicular axes to give different stability characteristics. For the orthotropic case, bending about the Y axis is referred to as “strong axis bending,” while bending about the X-axis is “weak axis bending” (Fig. 6(a)). While the circular straw can be bi-stable when bent about any axis, in Sec. 4.4, we will show that for orthotropic shapes, only the strong and weak axes have bi-stable states.

We investigate the bending stability as we increase d and change the cross section from isotropic to more orthotropic. Here, two different values of d are chosen to represent track shapes of different properties: d/r = 0.18 corresponds to a close-to-isotropic shape, while d/r = 0.75 is a cross section with clear orthotropy. The stability map consisting of 25 models with diverse combinations of t and θ₁ show how the bending stability changes with orthotropy. For the strong-axis bending of the track shape with d/r = 0.18, the transition boundary from mono-stability to bi-stability moves toward the left and top corner (Fig. 6(b), top). Compared to circular cross sections, more combinations of t and θ₁ acquire the second stable state by bending about the strong axis of track shapes. However, bending about the weak axis is bi-stable with less combinations of θ₁ and t (Fig. 6(b), the bottom plot). These trends are reinforced by increasing the orthotropy with d/r = 0.75 (Fig. 6(c)), where strong-axis bending offers more bi-stable configurations while weak-axis bending offers less. For instance, the model of t/R = 0.0175, θ₁ = 20 deg is bi-stable in the strong-axis bending process and mono-stable in weak-axis bending. In summary, when morphing circular cross sections to the track shape, strong-axis bending will possess a larger domain of bi-stable parameters t and θ₁, whereas the bi-stable domain of the weak-axis bending will decrease.

In addition to the orthotropic stabilities, the energy barriers E_b and the second stable state ϕ_b also change as the cross section becomes more orthotropic (Fig. 6(d)). We study the model of t/R = 0.0145, θ₁ = 25 deg which remains bi-stable in both bending directions. The energy barrier is the highest at d = 0 (a circular cross section) and gradually decreases with the distance d in both bending directions. Compared to the circular cross section, the strong-axis bending (d/r = 1.1) requires less energy input to trigger the second stable state. Energy input to activate bending bi-stability in the weak-axis direction is about 7% higher as compared to the strong-axis direction. As the cross section becomes more orthotropic (d/r rises from 0 to 1.1), the second stable state drops from 4.3 deg to 1.9 deg for the strong-axis bending, whereas the weak-axis bending increases only slightly from 4.3 deg to 4.45 deg. The decrease in the strong-axis stable state is intuitive. The partial inversion shortens the compressed side of the tube by a finite amount, and this shortening results in smaller rotations for deeper cross sections (increased d/r).

In Sec. 4.1, we showed that by changing the cross section into the track shape, the structure can achieve stable rotations with less energy input. The track shape brings orthotropy to the bending of the structure, but the curvature of the cross section stays non-negative everywhere. Therefore, to explore the influence of cross-sectional negative curvature, we study the dumbbell shape (Fig. 1(f)) in bending about both X and Y axes (Fig. 7(a)). The negative curvature region is configured by the parameter λ_φ. Starting from λ_φ → ∞ (the track shape), the shape does not change linearly with the λ_φ. When λ_φ = 2.5, the shape has only a small negatively curved region, but the stability behavior changes substantially with a large portion of the bi-stable domain flipping into mono-stable for both bending directions (Fig. 7(b)). The transition boundary of both bending axes moves toward the right and bottom corner as λ_φ further decreases, and the mono-stable parameter domain is increased when λ_φ = 1.5 (Fig. 7(c)). Thus, introducing negative curvature into the cross section prevents the structure from bi-stable bending in both directions. The structure with larger negative curvature region (smaller λ_φ) requires smaller shell thicknesses or larger slant angle θ₁ to be bi-stable. The stability trend for simple orthotropic shapes observed in Sec. 4.1 remains valid, where the strong axis bending has a larger bi-stable domain.

Energy barriers also increase with the negative curvature in the cross section, i.e., decrease with the controlling parameter λ_φ (Fig. 7(d)). For the strong-axis bending, the track shape of d/r = 0.75 only needs $E¯b=1.11$ to activate the second stable state, whereas the dumbbell shape with λ_φ = 1.5 requires $E¯b=1.65$⁠, showing a 50% increase. The weak-axis bending also needs more energy input for the same comparison. The higher energy requirement is due to increased bending in the negatively curved region as will be discussed in Sec. 4.3. While the energy requirement for bi-stable bending is higher, models with the dumbbell shape achieve the bent stable state at a similar rotational angle ϕ_b as the track shape. For the strong-axis bending, the stable rotation stays at the same level when the cross section transitions from the track shape to the dumbbell shape. Stable rotation of the weak-axis bending drops slightly from 4.4 deg to 4 deg.

The cross-sectional influence on the bending stability can be explained by examining the energy redistribution within the structure, as it undergoes a bi-stable transition. In Fig. 4(b), we showed that a corrugation with larger θ₁ can release much more stretching energy during the partial inversion, making the total energy to be bi-stable. Both the simulation and the linkage model suggest that the compression side dominates the energy behaviors. Figure 8 shows the energy redistributions around the compression side, from the peak point to the valley point of the stretching energy curve. The redistributions of stretching and bending energy are plotted for three different cases: the strong-axis bending of the track shape (A and D); the weak-axis bending of the track shape (B and E); and the strong-axis bending of the dumbbell shape (C and F). Among all three cases, only the strong-axis bending of the track shape is bi-stable. When comparing strong- to weak-axis bending (Figs. 8(a) and 8(b)), there is a much larger release of stretching energy in the strong-axis case (−1.84 mJ versus −0.87 mJ). In the strong-axis bending, material on the compression side has a high confinement within the curved region, which leads to the much higher energy change (Fig. 8(a)). In the weak-axis bending, there is less confinement, and stretching energy is only released near the curved region. The amount of bending energy redistribution does not change much between the strong and weak directions (1.68 mJ versus 1.43 mJ), and thus, the weak direction becomes mono-stable.

As shown in Fig. 7, negative curvature in the cross sections resists bi-stable bending. For cross sections with negative curvature, while the bars on the positively curved region have confinement and release energy, the bars around the concave area (middle of the cross section) will gain some stretching energy during the partial inversion (Fig. 8(c)). Thus, the release of stretching energy in the cross section with negative curvature gets slightly lower (−1.74 mJ) due to the gain near the concave area. The concave area also shows a higher increase in bending energy, as compared to the track shape (Figs. 8(f) versus 8(d)). Because there is less stretching energy released and more bending energy gained (2.3 mJ), negatively curved cross sections are more likely to be mono-stable.

For a specific bending axis, the directional valley point of the energy profile (Fig. 4(a)) is identified as the second stable configuration ϕ_b. Unlike the circular cross section, the two orthotropic shapes exhibit non-uniform energy profiles for different bending directions. The valley point (Fig. 4(a)) of a specific direction can be unconditionally stable only if the valley energy is the local minima among those of the neighboring directions, and such a direction is the so called stable bending direction. Therefore, stable bending directions can be identified by locating local minima of the energy function for bending directions and rotations ϕ. Models with the track shape and dumbbell shape are bent along all directions that are denoted by polar axes (Figs. 9(a), 9(b)), and bending axes are perpendicular to the associated directions. The cross sections with d/r = 0.9, θ₁ = 30 deg and t/R = 0.01 are selected so that all directional energy landscapes are bi-stable. The negative curvature controlling coefficient for the dumbbell shape is set to be λ_φ = 1.8.

For both shapes, strain energies for all bending directions are plotted in a polar coordinate system (Figs. 9(c), 9(d)). The distance from the centroid represents the rotation ϕ in degrees, and the magnitude of energy is distinguished by colors. Among all the bending directions, valley points of the strong-axis (0 deg or 180 deg) and the weak-axis (90 deg or 270 deg) possess the smallest strain energy in their neighborhood. Therefore, the structure is unconditionally stable in these four bent configurations. If the structure is bent about one of the non-primary directions (e.g., 45 deg), it will still pass over a peak load and will ultimately settle in one of the four bent stable states (Figs. 9(c), 9(d)).

The four stable configurations can also be found by plotting the energy magnitude for the different directions at specific rotational angles (Figs. 9(e), 9(f)). Distance from the centroid in Figs. 9(e), 9(f) represents the magnitude of energy, and different lines represent different levels of rotation. For rotational angles of ϕ = 1 deg, 1.5 deg which are before the peak point, the strong axis gives the largest strain energy for the rotation, corresponding to the highest bending rigidity. At valley points, the strain energies of the strong axis and of the weak axis are lower than those of their nearby axes, suggesting again that these two axes give the unconditionally stable states. The peak energy (energy barrier $E¯b$⁠) for all directions is plotted with a dashed line, and the overall shapes are close to circular. For the track shape, the shortest distance from points on the peak energy curve to the centroid is $E¯b=1.89$⁠, and the largest distance is $E¯b=2.0$⁠, only 5.8% different. For the dumbbell shape, the largest difference among all energy barriers is only 1.3%. The computed circularities of both peak-energy distribution curves are greater than 0.99 (1 for a circle), implying a high isotropy of the peak energies, despite the orthotropic bending behaviors. This means that roughly the same amount energy input is needed to activate bi-stable bending in all directions, regardless of the direction in which the system is loaded.

Orthotropic bending is possible with the track shape, and dumbbell shape cross sections, where bending along opposing directions returns identical responses. Therefore, a clover-like shape is introduced as a cross-sectional geometry to explore non-symmetric bending (Fig. 1(g)). To identify stable bending directions, a two-unit corrugation with the clover cross sections is examined by bending along all directions (Fig. 10(a)). The structure is defined by choosing θ₁ = 40 deg, t/R = 0.01 and d/r = 1.5, λ_φ = 1.5 to ensure all directional energies have bi-stable profiles. A tube with the clover shape cross section is expected to exhibit different behaviors when bending along the 0 deg or the 180 deg directions, while the behavior of bending along 90 deg and 270 deg directions should be identical.

Strain energies for all bending directions are plotted against rotational angles ϕ (Fig. 10(b)). The distance from the centroid represents the rotation ϕ in degrees, and the energy magnitude is differentiated by the different colors. The strain energy contours shows a rotational symmetry of order 3, inherited from the symmetry of the clover shape. Bending along 0 deg and 180 deg direction are both bi-stable, but the 0 deg-bending has a higher energy barrier and valley energy (Fig. 10(b)). Valley points of the above two bending directions are also local minima of the strain energy function, suggesting that these are unconditional stable states. Thus, the 0 deg and 180 deg directions are two stable bending directions, and rotational symmetry gives additional groups of stable directions: 120 deg, 300 deg and 60 deg, 240 deg. The other bending directions, for example 90 deg or 270 deg, do not have a local minima, and unrestricted bending in those directions would lead the structure to settle into one of the two nearby unconditionally stable states (e.g., 60 deg or 120 deg for the 90 deg bending).

The stabilities and symmetry of the bending directions can also be recognized by plotting the strain energies at specific deformation states (Fig. 10(c)). For specific rotations of ϕ = 1.2 deg,1.8 deg which are before the peak energy, directional strain energy grows with the bending depth. Bending along the stable directions of 60 deg, 180 deg, 300 deg, one has to overcome the smallest peak energy of $E¯b=9.4$⁠. The opposite directions have the largest peak energy of $E¯b=14.4$⁠, which is 53% larger. Comparing to the nearly isotropic energy barriers of the track and dumbbell variants, the clover shape shows a bigger difference among directional energy barriers. Strain energy of stable directions 60 deg, 180 deg, 300 deg slightly decrease to $E¯b=7.85$ (only a 16.5% decrease) at the valley, while these of their opposite directions drop by 30% to $E¯b=9.36$⁠. The valley energy curve in Fig. 10(c) also confirms that there are six stable directions, where the valley energies are local minima within their neighborhood.

Following the identification of stable bending directions, we explore the effects of the clover cross-section on the bending bi-stability with respect to geometric parameters: slant angle θ₁ and shell thickness t/R. We define the clover cross-section with a controlling coefficient λ_φ = 1.5, and we evaluated two values of the d/r ratio. We then explore the behavior for three directions of bending including: the 0 deg (Fig. 10(d)) and the 60 deg (Fig. 10(e)) directions which we know provide unconditionally stable states, and the 90 deg (Fig. 10(f)) direction which can have a bi-stable energy curve, but not an unconditional stable state. For the range of θ₁ ∈{20 deg, 25 deg, 30 deg, 35 deg, 40 deg} and thicknesses t/R ∈ {0.01, 0.0115, 0.013, 0.0145, 0.016}, all cases are bi-stable for bending in the 0 deg direction when d/r = 0.75. When the ratio d/r is increased to 1.5, there are a few less bi-stable bending cases for this direction (Fig. 10(d)). Bending in the direction of 180 deg, however, shows smaller domains of bi-stable parameters for both d/r ratios (Fig. 10(e)). These results show that, without switching the bending axis, some structures can be bi-stable in one direction of bending, and mono-stable if the bending direction is reversed (e.g. θ₁ = 25 deg and t/r = 0.13). Bending in the perpendicular direction (90 deg), which is unstable based on the directional energy analysis, is more sensitive to the distance ratio d/r (Fig. 10(f)). Among all twenty-five combinations of θ₁, t, almost all cases have a bi-stable energy curve with a smaller distance ratio d/r = 0.75. However, the number of bi-stable combinations drops significantly as the distance ratio increases to d/r = 1.5. For all three directions, increasing d/r ratios is a negative factor on the bending bi-stability, and the bi-stable domain shrinks with an increase in the d/r ratio. The result here is similar to that of the orthotropic shape with negative curvature. While the distance of the clover leafs increases, the shape with d/r = 1.5 has more negative curvature in the cross section which leads to less bi-stability.

This section describes physical models with the conventional circular and the track cross-sectional shapes made from paper sheets that serve as proof of concepts for some of the bending bi-stability behaviors (Fig. 11). Each of the specimens is created to have two unit cells to capture the interactions between cells (Fig. 2(c)). For each frustum, the planar projection is input into a CAD software and a paper sheet is cut using a laser cutter (Universal Laser VLS6.60). Each planar cut of the frustum has extra tabs for connecting them to the adjacent frustum to form a single unit cell. The same R/r ratio is used as in Sec. 2.1, where R = 31 mm and r = 27 mm are used for the paper models. The top frustum angle θ₁ is chosen as 30 deg, and d/r = 0.9 is selected for the track shape. The paper sheets are 0.25 mm (0.01 in) thick, and the normalized thickness is t/R = 0.008. The bar and hinge model is employed to investigate strain energy profiles.

The folding stiffness of the creases connecting frusta cannot be treated as zero, since the connection tabs have a finite value of crease stiffness. After calibration with respect to bending stabilities, L*/R = 3 is chosen to generate consistent predictions for all three bending cases in Fig. 11.

Here, three bending cases are studied: (1) circular cross-section; (2) the weak-axis bending of the track shape, and (3) the strong-axis bending the track shape. Starting from the initial state at which all unit cells are extended, these models are deformed by applying compression on one side of the system by hand (Fig. 11). After releasing the applied pressure, the model with circular cross sections can stay at the deformed state, implying a bi-stable bending behavior (Fig. 11(a)). The weak-axis bending of the track shape of d/r = 0.9 returns to the initial configuration and is mono-stable (Fig. 11(b)), whereas the strong-axis bending can maintain the bi-stable deformed configuration (Fig. 11(c)). The corresponding energy profiles for the three bending cases simulated with the bar and hinge model predict the stability consistently, using a L*/R = 3 to model folding stiffness (Fig. 11). These results confirm that for certain track shape geometries strong-axis bending can retain bending bi-stability, while the system becomes mono-stable for weak-axis bending.

Fabrication using paper models can introduce secondary effects that influence the bending stability and which are currently not explored with the bar and hinge simulation. Bending flat sheets to curved surfaces will introduce pre-stress that was previously found to be relevant to the bending stability [18]. The fabrication of the paper models may also introduce local defects, or introduce plasticity within the system. While our physical models demonstrate the effect of the geometry, we have not been able to investigate the influence of these secondary effects. We expect future research can give more insight to how such physical properties of the frusta effect the bending stability.

Thin-walled corrugated tubes made from open top frusta can allow for multi-stable global bending through partial inversions of their constituent cells. The capability to bend adaptively allows the tubes to conform into different configurations and to connect points in difficult to navigate environments. However, the bending bi-stability is not guaranteed for all corrugated tube designs. In this work, we investigate how the bending stability is affected by the cross-sectional geometry of the tubes and the different parameters that define the frusta.

We introduce and use a simplified bar and hinge model to capture the global behavior of different corrugated tubes. This model provides rapid and reasonably accurate predictions of the behavior with few convergence issues, which makes it well suited for the extensive parametric studies performed in this work. The bar and hinge model is calibrated and verified with respect to a more detailed, but much slower FE model. To understand and verify the most notable bending stability behaviors, we also explore representing the frusta using a simple four-bar linkage model where the sheet stretching and bending are represented intuitively. Finally, a physical paper model demonstrates how the bending stability of the corrugated tubes can be programmed by the cross-sectional geometry.

For tubes with a circular cross section, we find that a stable bending state can be achieved with higher conical angles of the top frustum, thinner shells, and a lower stiffness of the creases that connect separate frusta. All three measures increase the relative amount of stretching energy to bending and folding energy. Stretching energy is the only one that has a bi-stable profile, and is the driving component for achieving stable bent states of the tube.

Next, we investigate how the cross-sectional shape affects the bending stability. We first modify the circular cross section by moving two semicircles away to create a track shape, and then we create a dumbbell shape by squeezing the sides of the track shape. These doubly symmetric shapes can only achieve bent stable states in the two orthogonal directions. We show that bending bi-stability becomes more common for the strong axis bending as we increase the depth of the cross sections. When the track shapes are bent about the strong axis, the inverted region is curved and more confined so more stretching energy is released to cause a bistable bending. In contrast, for weak axis bending the inverted region is elongated with less stretching and reduced bi-stability. This anisotropic behavior of the track shape is confirmed with simple paper prototypes. Increasing the concavity of the cross section in the dumbbell shape leads to more mono-stable behavior for both directions of bending. The bending energy which increases monotonically becomes concentrated around the negatively curved region and counteracts the release of stretching energy to make the strain energy profile mono-stable. To further explore anisotropic bending, we evaluate a clover cross section that can have six stable bending directions with different stability behaviors when bending in opposing directions.

This paper gives an understanding to the geometry driven bending stability of thin-walled corrugated tubes and can inform future design of these systems to achieve desired multi-stable behaviors. We have characterized the stable bending angles which can be used to design reconfiguring multi-unit corrugated tubes that bend into desired complex shapes. The orthotropic cross sections have distinct directional stabilities, which allow for programmable bending with potentially two, three, four, or six stable bending directions that depend on the cross-sectional design. These systems can be used for tubes with predictable bending paths or robotic components that bend and re-orient in a desired fashion. The geometric design also influences the energy barriers for bi-stable bending, which can be used to design systems that only snap into a bent configuration after a specific force or moment is applied. Corrugated tubes with anisotropic cross-sections offer numerous areas for further research. For example, appropriate fabrication methods for such systems should be explored, and the strength of these tubes for orthogonal loads should be characterized. Given that the corrugated tubes can store fluids, it would also be interesting to investigate how internal pressure changes the stability behavior. The models and concepts introduced in this paper can be of use to this future research.

The authors acknowledge partial support from the National Science Foundation through Grant No. CMMI 1943723 and from the Office of Naval Research through grant N00014-18-1-2015 (# 12461330).

The authors report no conflict of interest.

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