Localized patterns in crushed conical shells

Buckling patterns in thin plates and shells arise due to the relative ease of bending a thin sheet compared to stretching it, and take a variety of forms, from global wrinkling to local folding and crumpling [1–3]. A particularly dramatic version of these patterns is seen in the buckling of axially loaded thin-walled cylinders; the reader can attest to this by standing on an empty aluminum can that will suddenly collapse into a crumpled, flattened disordered shape wherein the cylinder deforms on the largest scale possible, folding on itself [4]. This observation is in sharp contrast to theoretical studies which show that the energetically preferred patterns are a periodic array of buckled rhombi cells of typical size $(rt)^{1/2}$ (where r is the cylinder radius and t its thickness) composed of triangular planar surfaces known as the Yoshimura diamonds [5,6]. Indeed, decades of careful experiments confirm the results of the casual experiment alluded to earlier, i.e., producing ordered structures by compression of a cylinder is almost impossible. This is because the buckling instability of shells is well known to be subcritical, and thus very sensitive to imperfections [4,7,8]. Small imperfections associated with geometric defects thus grow rapidly on loading, and the entire structure collapses on large length scales with little reproducibility. However, it is possible to obtain the ordered Yoshimura patterns in experiments, a fact known to and exploited by artists and architects [9]. This requires one to prevent the growth of defects on large length scales, and can be done by introducing a slightly smaller rigid inner cylindrical mandrel inside the cylindrical shell before it is compressed [10] —this prevents the large wavelength buckling patterns from growing, so that the diamond-like ordered patterns arise naturally; the reader may have observed these patterns in tight-fitting trousers (or long soft boots) at the back of the leg, especially in the popliteal region. Given these observations that long wavelength subcritical instabilities are easy to excite in slowly loaded shells suggests that if a cylinder or cone is rapidly compressed to speed up the collapse, this would only amplify the disorder in the final crumpled state, in analogy with such systems as crystallization where slow annealing allows the material to solidify to an ordered, energetically favored configuration, while quick quenching drives the system into a disordered frustrated state.

Here we show that exactly the opposite is true —when a conical shell is smashed into a wall rapidly, it leads to the creation of an ordered 3D structure on small scales similar to inextensional fold patterns on cylinders and cones [4,11]. Our experiments use paper cones made by rolling flat paper sheets of thickness $t = \mathcal{O}(100\ \mu \text{m})$ and height $\mathcal{O}(20\ \text{cm})$ into cones of varying opening angles, with an overlap of the sheets such that each cone is two paper-layers thick. The shape of the cone is maintained by a double-sided tape along one of its generators at the edge of the sheet. The cones are mounted on a heavy circular base with the tip facing down so that the air is sealed within the cone, and dropped along a vertical rail onto a flat glass surface at velocities of $u = \mathcal{O}\ \text{(m/s)}$. A high-speed camera (Phantom V7.3) is used to image the cones from the side as they collapse upon impact with the surface at a frame rate of $\mathcal{O}(5\ \text{kHz})$.

As the tip of the paper cone impacts the rigid surface, a dynamical cascade of folds is initiated as seen in fig. 1(a) and Supplementary Movie SupplementaryMovie1.mp4. The cone deforms at the impact zone into localized rhombic patterns even as new rhombic indentations nucleate typically between two indentations of the previous generation. Although the cone deforms most strongly along a narrow circumferential rim in the neighborhood of the contact zone, we see that the cone is also weakly indented away from this zone. When the cone is unfolded by cutting it open along a generator and flattening it, we see that the process leads to an ordered tiled pattern —see fig. 1(b), (c). The rhombi composing the pattern may be characterized in terms of their azimuthal and radial length, $l_\theta$ and l_R, respectively, as shown in the inset to fig. 1(c). l_R initially increases with the radial coordinate R, but quickly saturates at a constant value as shown in fig. 2(a). In contrast, $l_\theta$ continues to increase well beyond the point at which l_R saturates, leading to more elongated cells farther from the cone tip.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The form of the patterns depends on a number of parameters such as the cone opening angle, paper thickness, and impact velocity, and so we systematically vary them. Although the ordered patterns only form under dynamic conditions, small variations in the velocity (by a factor of three) do not change the observed pattern, as seen in fig. 2(b). Changing the cone opening angle, $\alpha=\tan^{-1}(r/R)$, where r is the local cone radius indicated in fig. 2(c) changes the form of the rhombi. For larger α, $l_\theta$ grows more rapidly with R, as shown for four different opening angles in fig. 2(c). This suggests a dependence of $l_\theta$ on r, as shown by the data collapse onto one master curve in the inset to fig. 2(c), indicating that the shape of the pattern is determined locally by the instantaneous value of r in contact with the surface. In contrast to the strong dependence of $l_\theta$ on the cone geometry, we do not measure any dependence of l_R on α.

To characterize the role of sheet thickness on the pattern, we varied it from the typical value we used, 0.1 mm, to a range spanning 0.075–0.28 mm. We find that the effect of thickness on the length scales of the rhombi pattern was very small. This observation strengthens the claim that the geometry of the crushed pattern is governed by the dynamics during the crushing process, rather than by minimization of elastic energy. However, the paper thickness does affect the region of the cone where a Yoshimura-like pattern emerges. At the very tip of the cone, where the rhombi sizes are expected to be very small, the paper creases instead of buckling into a rhomboidal pattern over a distance that is strongly affected by thickness, from few millimeters for 0.075 mm thick paper to about 5 cm for 0.28 mm thick paper.

While small variations in the velocity of impact do not change the observed folded patterns, when the cone is compressed quasi-statically the patterns disappear and the cone crumples in a disorderly fashion with no localization (see Supplementary Movie SupplementaryMovie2.mp4). This difference suggests a role for dynamic effects associated with inertia in the shell and a potential role for air pressure in stabilizing the long wavelength modes excited by quasi-static compression. The relative role of inertial to elastic forces (per unit area) is determined by the ratio ${\sim}\frac{\rho_{paper} u }{\tau}/\frac{B }{L^3}$ [8]. Substituting parameter values ($t=100\ \mu \text{m}$, $\rho_{paper}=800\ \text{kg/m}^{3}$, $u=1\ \text{m/s}$, $B=E t^3$ with Young's modulus of paper $E=2\ \text{GPa}$ and the typical length and time scales for buckling rhombi $L=1\ \text{cm}$ and $\tau=0.01\ \text{s}$, respectively), we find that the ratio ${\sim}10^{-2}$, i.e., inertial effects are negligible. This leaves only the air pressure as a pattern stabilizer, wherein it plays the role of a dynamic mandrel to suppress the growth of long wavelength defects in favor of the short wavelength ordered patterns that form at the zone of impact. To test this experimentally, we mount the cones on bases through which four circular holes of varying diameter, d, are cut, as shown schematically in fig. 3(c). This allows air to escape the cones and reduce the rate of buildup of the dynamic pressure so that the compact patterning associated with tip collapse becomes less ordered —fig. 3(a).

Zoom In Zoom Out Reset image size

To quantify this change in pattern evolution as a function of the hole size (see Supplementary Movies SupplementaryMovie1.mp4 and SupplementaryMovie3.mp4), we measured the size of the indentation zone h shown in fig. 3(b) as a function of the axial distance H from the tip of the cone for different hole sizes, where h and H are extracted from real time imaging of the crushing process closely and correlate with $\ell_R$ and R, respectively. For small holes, the air pressure inside the cone is high, and the new indentations remain localized near the interface of the cone and the surface. For larger holes, the air pressure inside the cone is small, and once indentations are created they grow rapidly toward the base of the cone. In addition to the increase in size the indentations reach, the release of air pressure also increases the chance for a dent to nucleate away from the surface, causing the cone to buckle catastrophically (see Supplementary Movie SupplementaryMovie4.mp4). In the absence of holes, the increase in the internal pressure $\Delta P$ scales as the ratio between the uncompressed volume of the cone and the truncated compressed cone, of half the volume, and yields $\Delta P \approx 1\ \text{atm} \approx 100\ \text{kPa}$. Direct measurement during compression confirms this buildup of pressure within the cone; we observe a peak pressure difference of between 16 and 27 kPa. When the holes are large enough so that as the air escapes viscous and compressible effects are negligible, $\Delta P \approx \rho_{air} u_{air}^2$, where $\rho_{air}$ and u_air are the density and velocity of the escaping air, respectively. The velocity u_air is estimated by equating the flux of air escaping the cone to the rate of the cone's volume change. For a cone impacting the surface at a velocity of 1 m/s, this calculation yields $\Delta P \approx 100\ \text{Pa}$, three orders of magnitude lower than in the pressurized case, confirming that the dominant stabilizing effect that causes localized folding at the impact zone is due to pressurized air in the cone.

We validated the buildup of pressure by measuring the internal gauge pressure as the cone compacts with an ELVEFLOW MFP Inline Pressure Sensor FlowPlus embedded in the top circular plate sealing the cone. The cones were dropped from a height of 1.2 m and pressure data recorded at a rate of 500 Hz. The results are plotted in fig. 4; the interior of the cone attains a peak internal pressure of between 16 and 27 kPa. Leaks at the interface between the top plate and the cone's base, porosity of the paper, and leaks due to delamination of the cone along its lateral seam all contribute to the fall in pressure at the end of the experiments, and the discrepancy between the measured peak pressure and predicted theoretical peak pressure of 100 kPa.

Zoom In Zoom Out Reset image size

To understand our results quantitatively, we recall that an energy-based linear stability analysis predicts that the typical length scale of the buckling pattern would scale as $l\sim\sqrt{r t}$. However, the buckling patterns observed do not follow this scaling relation because the pattern evolves dynamically and cannot be described purely in terms of energetics. Therefore, we numerically simulated compression of conical shells using a thin-shell theory based on the Kirchhoff-Love assumptions that the shell geometry and mechanics can be characterized entirely in terms of the geometry of its midsurface. Furthermore, we consider only geometric nonlinearities, and limit ourselves to elastic buckling since this always precedes plastic deformations and determines the basic patterns.

The shell can be parameterized entirely by the geometry of a midsurface S; more specifically, let A be an annulus in the plane and $r:A \rightarrow {R}^3$ the coordinate function embedding S (illustrated schematically in fig. 5(a)). The thin shell of thickness t is then parameterized by

$$\begin{equation*} \tilde{r}: A\times[-t/2,t/2] \rightarrow {R}^3, (x,y,z)\mapsto r(x,y) + z\hat{n}(x,y), \end{equation*}$$

where $\hat{n} = (r_x\times r_y)/\|r_x \times r_y\|$ is the surface normal of S. For the large deformations of a thin conical shell, geometric nonlinearities dominate material nonlinearities, and so we use a simple linear St. Venant-Kirchhoff material model. Denoting by $\bar{S}$ the embedding of the cone at rest, before compression, and letting $\bar{r}$ be its corresponding embedding, the Steiner expansion formula [12] then gives the elastic energy in terms of geometric quantities on the midsurface: [13,14]

$$\begin{equation} E = \int_S\int_{-t/2}^{t/2}2k_1 \mathrm{Tr}[\epsilon(\eta)^2] + k_2 \mathrm{Tr}^2[\epsilon(\eta)]\,\textrm{d}\eta\,\textrm{d}\bar{A} \end{equation} \tag{ 1 }$$

where $\epsilon(\eta) = \bar{a}^{-1}(a-\bar{a})+2\eta\bar{a}^{-1} b + \eta^2 \bar{a}^{-1}c$, and a, b, c are the first, second, and third fundamental forms of S that can be expressed in terms of r by $a = \mathrm{d}r^T\mathrm{d}r$, $b = \mathrm{d}r^T\mathrm{d}\hat{n}$, and $c=\mathrm{d}\hat{n}^T\mathrm{d}\hat{n}$. Here $\mathrm{d}\bar{A}$ is the rest surface area element, and is discretized using a piecewise-constant metric over S, $\bar{a}$ is the rest metric of the midsurface (the first fundamental form $\mathrm{d}\bar{r}^T\mathrm{d}\bar{r}$ of $\bar{S}$), and k₁, k₂ are Lamé parameters.

Zoom In Zoom Out Reset image size

Since we expect the in-plain strain to remain small during compression, we discretize S using piecewise-linear triangular elements $\mathcal{T}$; following Weischedel [15], we discretize the fundamental forms and write

$$\begin{equation*} E \approx \sum_{T \in \mathcal{T}} A(T) \int_{-t/2}^{t/2}2k_1 \mathrm{Tr}[\epsilon(\eta)^2] + k_2 \mathrm{Tr}^2[\epsilon(\eta)]\,\mathrm{d}\eta, \end{equation*}$$

where A(T) is the area of triangle T. To determine expressions for the fundamental forms we start by assuming that q_i, for $i\in\{0,1,2\}$, denote the positions of the triangle vertices, $\bar{q}_i$ the corresponding vertex positions on the undeformed cone S, and p_i their preimage in the parameter domain A. We can choose coordinates on A so that $p_0 = (0,0), p_1 = (1,0), p_2 = (0,1)$ so that the first fundamental form

$$\begin{equation*} a = \left[\begin{array}{c@{\quad}c} \|q_1-q_0\| & (q_1-q_0)\cdot(q_2-q_0)\\ (q_1-q_0)\cdot(q_2-q_0) & \|q_2-q_0\| \end{array}\right] \end{equation*}$$

with $\bar{a}$ defined similarly, with $\bar q_i$ taking the role of the q_i. Discretizeing the second fundamental form by way of midedge normals: we assign a unit vector $n_{e_i}$ to the midpoint of edge e_i opposite q_i, where $n_{e_i}$ is equal to the face normal, if e_i is on the top or bottom boundary of the cone, or the normalized average of the two incident triangles, if e_i is an interior edge. Since the segments connecting triangle edge midpoints are parallel to and half the length of the triangle sides, we can write

$$\begin{equation*} b = 2\left[\begin{array}{c@{\quad}c} (n_{e_0}-n_{e_1})\cdot(q_1-q_0) & (n_{e_0}-n_{e_1})\cdot(q_2-q_0)\\ (n_{e_0}-n_{e_2})\cdot(q_1-q_0) & (n_{e_0}-n_{e_2})\cdot(q_2-q_0) \end{array}\right]. \end{equation*}$$

Note that while $a, \bar{a}$ and b depend on the choice of coordinates on A, the energy E does not. The above discretization can be interpreted as the discrete shape operator [16] in the case where the edge directors are fixed to equal the average face normals. Finally the third fundamental form c is computed from a and b using the identitity $c = 2Hb - Ka$, where $H = \frac{1}{2}{\mathrm{Tr}}(a^{-1}b)$ and $K = \mathrm{Det}(a^{-1}b)$ are approximations of the mean and Gaussian curvature of S.

Since the cone is subject to body forces due to the dynamic buildup of pressure, we add an external force to the vertices proportional to the pressure difference $\Delta P$ across the cone; the force applied to vertex i is $F_{\mathrm{pressure}} = A_i \Delta P\,n_i$, where the vertex barycentric dual area A_i is one-third the sum of the faces incident to vertex i. For the discretization of the surface normal n_i at the vertex, we chose the mean curvature normal, based on the relationship $\Delta r = 2Hn$, where Laplace-Beltrami operator Δ is applied component-wise. This relationship is discretized to yield $n_i = (Lp)_i/\|(Lp)_i\|$, where L is the cotangent Laplacian [17] and p is the vector of vertex positions.

Crushing of the cone was modeled by simulating a plane of given mass density, initial height H, and initial downward velocity, and adding a unilateral inequality contact constraints between the plane and every vertex of the cone mesh. Both gravity and contact forces act on the plane, so that it initially accelerates, then decelerates as pressure builds up in the cone. We neglected self-contact as it does not occur until after the creases of the rhombus pattern have formed. In a typical simulation, the entire cone is composed of roughly 160000 elements, and the evolution of the cone is integrated in time using a velocity Verlet algorithm with a time step size of $2\times 10^{-7}$ seconds. To damp the transient elastic waves that propagate through the cone during compression, we imposed a Rayleigh damping model with a damping coefficient of 0.001. We use 160000 triangular elements to discretize the cone and write the elastic energy in terms of the discrete versions of the first and second fundamental forms [15,16,18]. We apply an external force to the vertices proportional to the pressure difference $\Delta P$ across the cone, and then solve the discrete equations of motion using a velocity Verlet scheme using a time step of $2\times 10^{-7}$ seconds. Our simulations allow us to vary α and the pressure inside the cone, and measure how these parameters affect the pattern as well as the elastic energy stored in the deformed cone.

For the case of a sealed cone, we approximate the pressure within the cone as being inversely proportional to the volume ratio between the crushed and uncrushed cone. When simulating the crushing of cones with holes allowing for the release of pressure, we estimate the pressure as follows. When the holes in the cone's base are large, air escapes quickly through the opening and the pressure within the cone is expected to be negligible relative to that within a sealed cone (no holes in the base). We assume that in the large-hole regime, viscous and compressible effects are negligible, and $\Delta P \approx \rho_{air} u_{air}^2$, where $\rho_{air}$ and u_air are the density and velocity of the escaping air, respectively. Incompressibility allows us to compute u_air based on the geometry of the cone and its compression velocity u_cone: the volume air escaping through the holes during an increment of time dt is $n_{holes} \pi r_{holes}^2 u_{air}\mathrm{d}t$, where n_holes and r_holes are the number of holes in the base of the cone, and their radii, respectively. The change in volume due to compression of the cone is $\pi r(H)^2 u_{cone}\mathrm{d}t$, where r(H) is the radius of the portion of the truncated cone in contact with the ground, $r(H) = r_0 H/H_0$. Equating the volume loss due to compression and the flux of air through the holes yields, for the case of four holes, $\Delta P \approx \rho_{air} \frac{u_{cone}^2 r_0^4 H^4}{16 r_{holes}^4 H_0^4}$. This formula is used to estimate a peak pressure difference of $\approx100\ \text{Pa}$ for a cone with holes of diameter 17 mm, and is used in the numerical simulations to evaluate the internal pressure of the cone at every time step.

In fig. 5(b) we show that our simulations capture the rhombic patterns seen in our experiments (see Supplementary Movie SupplementaryMovie5.mov), despite the lack of plasticity in the simulated sheets [19–21]. In fig. 5(c), we see that the dependence of l_R and $l_\theta$ on R observed in the simulation are in agreement with the experimental data. However, the boundaries between neighboring rhombi appear to be less sharp and well defined in simulations than in experiments, which we attribute to the simple elastic constitutive law used in the simulations.

For the simulations using pressurized cones, the P – V work done by the air on the walls of the cone becomes significantly larger than the elastic deformation energy stored in the sheets, as can be seen in fig. 5(d). When we simulate compression of a cone where air pressure is released via four holes of diameter 17 mm, the energy required to compress the air drops by three orders of magnitude, well below the elastic energy for all values of compression, indicating that air pressure no longer plays a role in the deformation of the cone.

In fig. 6 we demonstrate the effect of pressure on the buckling process in simulation. For high pressure, the deformations localize to the interface with the surface (fig. 6(a)), while for low pressure, deformation of the shell extends far away from the interface (fig. 6(b)). The most notable difference between the simulations and experiments is that the deformations observed in simulations with low pressure, while not localized, still appear as ordered layers of rhombi extending away from the surface, as opposed to a smaller number of large buckles. We attribute this to the ability of simulation to find better local minima in the energy landscape of the deformed shells, and believe that in the noisy experimental environment the lack of localization results in the large disordered buckling patterns we observe. Despite this change, we do not observe an increase in the length scales of the buckled rhombi or a transition to disorder as the pressure is released. We attribute this to the lack of large material and geometrical inhomogeneities which are known to be important for the nucleation and propagation of defects [4,7,8].

Zoom In Zoom Out Reset image size

Overall, our study shows how it is possible to suppress subcritical instabilities in a rapidly compressed conical shell by taking advantage of the role of air pressure that works like a dynamic mandrel, and prevents the growth of long wavelength instabilities. This in turn forces the localization of pattern to small scales in the vicinity of the impact zone from where it builds up hierarchically. More generally, this is suggestive of a general principle by which a dynamical process may be used to prevent the growth of (slow) long wavelength instabilities and thus lead to localized patterns and order in systems susceptible to subcritical instabilities.