Foldable Cones as a Framework for Nonrigid Origami

TL;DR

This paper investigates the elastic behavior of foldable cones (f-cones) in origami, incorporating nonlinear effects, crease mechanics, and testing inextensibility, offering a new numerical modeling approach for realistic origami structures.

Contribution

It introduces a comprehensive nonlinear elastic model for f-cones, including crease mechanics and deformation effects, advancing beyond previous linear models.

Findings

Nonlinear corrections significantly affect f-cone behavior.

The phase field-like model effectively simulates realistic origami deformations.

Inextensibility hypothesis is validated through numerical testing.

Abstract

The study of origami-based mechanical metamaterials usually focuses on the kinematics of deployable structures made of an assembly of rigid flat plates connected by hinges. When the elastic response of each panel is taken into account, novel behaviors take place, as in the case of foldable cones (f-cones): circular sheets decorated by radial creases around which they can fold. These structures exhibit bistability, in the sense that they can snap-through from one metastable configuration to another. In this work, we study the elastic behavior of isometric f-cones for any deflection and crease mechanics, which introduce nonlinear corrections to a linear model studied previously. Furthermore, we test the inextensibility hypothesis by means of a continuous numerical model that includes both the extended nature of the creases, stretching and bending deformations of the panels. The results show that this phase field-like model could become an efficient numerical tool for the study of realistic origami structures.

Figures

Equations

$${\mathbf{u}}^{\prime }$$

$${\mathbf{t}}^{\prime }$$

$${\mathbf{n}}^{\prime }$$

$$\sum^{n}_{i=1}\Delta\varphi_{i}=\left\{\begin{array}[]{ll}\pm 2\pi&\text{if }\Gamma\text{ encloses the }z\text{-axis}\\ 0&\text{if not}\end{array}\right.,$$

$$F_n[\mathbf{u},\mathbf{t}]$$

$${\kappa }_i^{\prime \prime }+\left(1+c_i\right){\kappa }_i+\frac{{\kappa }_i^3}{2}=0,$$

$$a\displaystyle\sum^{n}_{i=1}\left[-(\mathbf{f}_{i}\cdot\delta\mathbf{u}_{i})+(\kappa_{i}\mathbf{n}_{i}\cdot\delta\mathbf{t}_{i})\right]\Big{|}^{s_{i+1}}_{s_{i}}+\displaystyle\sum^{n}_{i=1}\delta g_{i}=0,$$

$${\kappa }_i^{+}={\kappa }_i^{-},$$

$${\mathbf{f}}_i^{+}={\mathbf{f}}_i^{-}.$$

$${\kappa }_i^{\prime +}$$

$$\kappa_{i}(s)=\left\{\begin{aligned} &\kappa_{0i}\,\text{cn}\left(\frac{\kappa_{0i}}{2\sqrt{m^{\mbox{\tiny R}}_{i}}}s-S_{0i}\Bigg{|}m^{\mbox{\tiny R}}_{i}\right),&&\text{if}\ J^{2}>c^{2},\\ &\kappa_{0i}\,\text{dn}\left(\frac{\kappa_{0i}}{2}s-S_{0i}\Bigg{|}m^{\mbox{\tiny S}}_{i}\right),&&\text{if}\ J^{2}<c^{2},\end{aligned}\right.$$

$$m_i^{\text{mbox}R}=\frac{1}{m_i^{\text{mbox}S}}=\frac{{\kappa }_{0i}^2}{2{\kappa }_{0i}^2+4(1+c)}.$$

$$g_i=L\left({\sigma }_i{\mathbf{t}}_i^{-}\cdot {\mathbf{t}}_i^{+}+{\tau }_i\left({\mathbf{t}}_i^{-}\times {\mathbf{t}}_i^{+}\right)\cdot {\mathbf{u}}_i^{+}\right),$$

$${\sigma }_i=k_i\cos {\psi }_i^0,{\tau }_i=k_i\sin {\psi }_i^0,$$

$${\kappa }_i^{+}=\frac{1}{a}\frac{d{g}_i}{d{\psi }_i}.$$

$$\delta F_n$$

$${\mathbf{f}}_i^{\prime }=({\kappa }_i^2-{\lambda }_i){\mathbf{u}}_i-a{\kappa }_i{\mathbf{n}}_i,$$

$${\mathbf{f}}_i={\kappa }^{\prime }{\mathbf{n}}_i+2{\kappa }_i^2{\mathbf{t}}_i-{\Lambda }_i{\mathbf{t}}_i.$$

$${\mathbf{f}}_i^{\prime }\cdot {\mathbf{t}}_i$$

$${\mathbf{f}}_i\cdot {\mathbf{u}}_i$$

$${\mathbf{f}}_i={\kappa }_i^{\prime }{\mathbf{n}}_i-\left(\frac{1}{2}{\kappa }_i^2+c_i\right){\mathbf{t}}_i.$$

$${\mathbf{f}}_i^{\prime }\cdot {\mathbf{u}}_i=\frac{1}{2}{\kappa }_i^2+c_i.$$

$${\mathbf{f}}_i^{\prime }=\left(\frac{{\kappa }_i^2}{2}+c_i\right){\mathbf{u}}_i-{\kappa }_i{\mathbf{n}}_i.$$

$${\mathbf{f}}_i^{\prime }=\left[{\kappa }_i^{\prime \prime }+{\kappa }_i\left(\frac{{\kappa }_i^2}{2}+c_i\right)\right]{\mathbf{n}}_i+\left(\frac{{\kappa }_i^2}{2}+c_i\right){\mathbf{u}}_i.$$

$$\mathbf{t}$$

$$\mathbf{n}$$

$$\mathbf{u}$$

$${\mathbf{n}}_{\phi }$$

$$\delta \mathbf{u}$$

$$\delta \mathbf{t}$$

$$\mathbf{n}\cdot \delta \mathbf{u}$$

$$\mathbf{t}\cdot \delta \mathbf{u}$$

$$\mathbf{n}\cdot \delta \mathbf{t}$$

$$\left(\begin{array}{c}{\mathbf{t}}_i^{+}\\ {\mathbf{n}}_i^{+}\end{array}\right)=\left(\begin{array}{cc}-\cos {\psi }_i& \sin {\psi }_i\\ -\sin {\psi }_i& -\cos {\psi }_i\end{array}\right)\left(\begin{array}{c}{\mathbf{t}}_i^{-}\\ {\mathbf{n}}_i^{-}\end{array}\right).$$

Full text

The basic premise of origami, the ancient Japanese art of paper folding, is to obtain a complex 3-dimensional structure starting from a 2-dimensional sheet to which a network of creases is imprinted. Despite the simplicity of this idea, in recent years, the field of mechanical metamaterials has sought inspiration from origami Schenk and Guest (2013); Wei et al. (2013) in the search of smart-materials with a vast range of functionality such as deployability of large membranes Miura (1985), shape changing structures Dias et al. (2012); Dudte et al. (2016), and tunable mechanical and thermal properties Silverberg et al. (2014); Waitukaitis et al. (2015); Boatti et al. (2017); Pratapa et al. (2019), just to name a few. In practice, many of these applications are constrained to situations which origami structures are made from assemblies of flat rigid plates connected by hingelike creases. In such situations, the geometrically accessible configurations are fully determined by the crease network, while the structural response is a result of the crease network and the crease mechanics Hanna et al. (2014); Brunck et al. (2016); Chen et al. (2015). By contrast, when the elastic response of the plates (mainly bending) is taken into account, a variety of new behaviors may emerge. In this case, the elastic response of the structure is determined by the competition between the flexural stiffness of the panels $B$ and the torsional rigidity of the creases $k$. The length $L^{*}\equiv B/k$, called origami length, determines wether the deformation of a non-rigid origami is bending or crease dominated Lechenault et al. (2014). If $l$ is the typical size of the facets, when $l\ll L^{*}$, the deformation is governed by the change on the folding angles, while if $l\gg L^{*}$, the deformation is governed by the bending of the panels.

However, suitable analytical models capturing the elastic regime of non-rigid origami still remains for the most part unexplored. Foldable cones Lechenault and Adda-Bedia (2015), or $f$-cones, are the simplest single-vertex non-rigid origami in which the elasticity of the plates is relevant. $f$-cones are elastic sheets decorated by straight creases meeting at a single vertex around which they are folded. As a first approximation, these sheets are assumed to be inextensible. This results in a family of various umbrella-like motifs whose equilibrium shapes depend on the crease pattern imprinted in the flat configuration of the sheet and the mechanical response of the creases. Regardless of the initial crease pattern, these structures exhibit bistability in the sense that they can mechanically snap through from one metastable configuration to another of higher elastic energy.

The $f$-cone belongs to a larger family of singularities emerging on sheets subjected to isometrical deformations Huffman (1976); Seffen (2016); Guven et al. (2011). In many situations, the elastic energy in a thin elastic sheet can localize in a single point leading to conical dislocation. The most fundamental example is the so-called $d$-cone Ben Amar and Pomeau (1997); Cerda and Mahadevan (1998), a conical singularity observed when crumpling an elastic thin sheet. The bending energy of the defect diverges logarithmically as one gets closer to the vertex. In a more realistic situation, these divergencies are regularized if the inextensibility constraint is relaxed, thus leading to stretching and plastic deformations close to the vertex Cerda et al. (1999).

The bistable behavior of $f$-cones was investigated in Lechenault and Adda-Bedia (2015) and a model was proposed to describe their equilibrium shapes in the limit of small deflection and infinitely stiff creases. We will refer their model as the linear model for $f$-cones, as it relies on the approximation of small deflections which allows to write the curvature of the surface as a linear function of the vertical component of the displacement. This system has motivated the study of other similar problems such as the bistable behavior of creased strips Walker and Seffen (2018). In this last work, the authors proposed a discrete model based on the Gauss map of several creases meeting at the vertex. In the limit of infinite creases, the linear version of an $f$-cone with two creases is recovered. The discrete model that is based on the Gauss map has limitations, as it only predicts the final shape of real sheets well for small deflections, while important discrepancies with experiments are observed for large deflections. Although, these discrepancies may be attributed to the existence of stretching in real sheets, which in turn invalidates the inextensibility hypothesis, the inherent non-linear nature of the system may also have a significant contribution to interpret the experimental observations. In the present work we propose an alternative model for $f$-cones that encompasses the full geometric non-linear contributions, thus capturing any deflections—this model describes the equilibrium shapes as function of the folding (dihedral) angle of the creases. Also, the effect of crease mechanics with hingelike behaviors is incorporated into the model. Then, we corroborate the predictions of the model with the aid of FEM simulations. From the numerical model we are able to quantify the stretching on the system in order to test the validity of the inextensibility hypothesis during the entire indentation process.

The manuscript is organized as follows. In section II the system under study and its geometry are presented in detail. Then, in section III, we present our elastic model for $f$-cones and the main results. The results presented here complement the predictions of the linear version of the model Lechenault and Adda-Bedia (2015). Subsequently, in section IV, a numerical model that simulate an $f$-cone of 4 creases is proposed to study the snapping process in a finite element analysis. Then, in section V the results of the numerical study are compared with the theory. The details of the analytical calculations can be found in the Appendix.

I Kinematical description of nonrigid single vertex origami

Foldable cones, or $f$-cones, are made from a circular elastic sheet decorated by one or more straight radial creases meeting at a single vertex Lechenault and Adda-Bedia (2015). These surfaces resemble those of $d$-cones Ben Amar and Pomeau (1997), except that they can fold around the creases. When the elastic sheet is inextensible, the only possible equilibrium shapes are developable surfaces and, in this particular geometry, developable cones. This implies that the deformed shape can always be isometrically mapped to the initial flat state. The equilibrium shape of the cone will be developable anywhere except at the tip of the cone and the creases, where the curvature is not defined. In this section, we first introduce the parametrization of a general conical shape, and then we describe in detail the geometry of an $f$-cone.

I.1 Geometry of developable cones

The most general parametrization of a conical shape is given by $\mathbf{r}(r,s)=r\mathbf{u}(s)$, where $r$ is the distance to the tip, $\mathbf{u}(s)$ is a unit vector and $s\in\left[0,2\pi\right]$ is the arc-length of the curve $\Gamma:s\rightarrow\leavevmode\nobreak\ \mathbf{u}(s)$ on the unit sphere. The tangent vectors adapted to the surface of the cone are $\mathbf{u}$ and $\mathbf{t}=\mathbf{u}^{\prime}$, where the prime denotes derivative with respect to $s$. As $s$ is the arc-length of the curve, the tangent vector $\mathbf{t}$ is a unit vector. Note that $\mathbf{u}\cdot\mathbf{t}=0$ and that the normal of the surface is given by $\mathbf{n}=\mathbf{u}\times\mathbf{t}$. Therefore, the triad $\left\{\mathbf{u},\mathbf{t},\mathbf{n}\right\}$ forms a right-handed basis that satisfies the following equations Cerda et al. (2004); Cerda and Mahadevan (2005)

[TABLE]

where $\kappa(s)=\mathbf{t}(s)\cdot\mathbf{n}^{\prime}(s)$. The metric tensor of the conical surface is given by $g_{ab}=\partial_{a}\mathbf{r}\cdot\partial_{b}\mathbf{r}$, where the indices $a,b=r,s$. Hence, for conical geometries, the metric components are $g_{rr}=1$, $g_{rs}=0$ and $g_{ss}=r^{2}$. The extrinsic curvature tensor is defined as $K_{ab}=\partial_{a}\mathbf{r}\cdot\partial_{b}\mathbf{n}$ and its single non-vanishing component is $K_{ss}=r\kappa$. Therefore, the surface curvature is $K=g^{ab}K_{ab}=\kappa/r$. Once $\kappa(s)$ is known, the final shape of the cone can be reconstructed by integrating Eqs. (1).

I.2 Geometry of foldable cones

We consider a $f$-cone of $n$ creases made from a flat circular sheet of radius $R$ which is parceled out in $n$ circular sectors (panels) delimited by the creases (see Fig. 1(a)). A hole of radius $r_{0}\ll R$ is cut out at the center in order to avoid a divergence in the elastic energy. Let $\alpha_{i}$ denote the sector angle of the $i^{th}$ panel in the flat configuration, with $\sum^{n}_{i=1}\alpha_{i}=\leavevmode\nobreak\ 2\pi$. The value of the arc-length at the $i^{th}$ crease is denoted by $s_{i}$, so that, through the inextensibility condition, $\alpha_{i}=s_{i+1}-s_{i}$ in the deformed configuration. Hereinafter, for any scalar or vector field of the form $b_{i}(s)$, the subscript $i$ specifies that the domain of the function corresponds to the $i$-th sector, where the periodic convention $b_{i\pm n}\equiv b_{i}$ is assumed. Moreover, we introduce the following notation: $b^{-}_{i}\equiv b_{i-1}(s_{i})$ and $b^{+}_{i}\equiv b_{i}(s_{i})$.

In the deformed configuration, each crease has a folding angle (dihedral angle) $\psi_{i}$, which we call it mountain if $\psi_{i}\in[0,\pi]$ or a valley if $\psi_{i}\in[\pi,2\pi]$. For a mountain and a valley creases, Fig. 2 shows how the folding angle is defined by $\mathbf{t}^{+}_{i}$ and $\mathbf{t}^{-}_{i}$ in the plane perpendicular to the crease. In terms of spherical coordinates, the final shape is given by a polar angle $\theta(s)$ and an azimuthal angle $\varphi(s)$ which are functions of the arc-length [Fig. 1(b)]. Each angular sector will span an azimuthal angle $\Delta\varphi_{i}\equiv\varphi(s_{i+1})-\varphi(s_{i})$. The closure condition can be written as

[TABLE]

where $\pm$ indicates that $\varphi_{i}$ could increase clockwise or counter-clockwise in the $x$-$y$ plane. Notice that in principle $\Delta\varphi_{i}$ and $\alpha_{i}$ are not necessarily equal in the deformed configuration. However, they coincide in certain symmetrical cases: $f$-cones with an arbitrary number of evenly distributed mountain creases or an even number of evenly distributed alternating mountain-valley creases where all the creases are identical. In such cases, we shall say in the following sections the $f$-cone is symmetrical.

II Elastic theory of foldable cones

Our model is based on a generalization of the functional introduced in Ref. Guven and Müller (2008). The total energy of an $f$-cones with $n$ creases is the sum of the elastic energy over all the panels plus the mechanical energy stored in the creases. Thus, the principle of virtual work is equivalent to minimizing the following functional

[TABLE]

Here, $a=B\ln\left(R/r_{0}\right)$, where $B$ is the flexural stiffness (bending modulus) of the sheet. The first term inside the brackets accounts for the bending energy of the facets, where $\mathbf{u}_{i}(s)\cdot(\mathbf{t}_{i}(s)\times\mathbf{t}^{\prime}_{i}(s))=\kappa_{i}(s)$ is the dimensionless curvature of the $i$-th panel. The above augmented energy functional contains $3n$ local Lagrange multipliers, namely $\lambda_{i}(s)$, $\Lambda_{i}(s)$ and $\mathbf{f}_{i}(s)$, which correspond to the following kinematical constraints, respectively: $\lambda_{i}(s)$ enforces $\mathbf{u}_{i}$ to be a unit vector, thus constraining the final trajectories to the unit sphere; $\Lambda_{i}(s)$ enforces the parameter $s$ to be the arc-length of the curve $\Gamma$; and finally, $\mathbf{f}(s)$ is a force (normalized by $a$) that anchors the tangent vector to the embedding. The functions $g_{i}$, which depend on the frame vectors at both sides of the crease, account for the elastic energy stored in the $i$-th crease. For simplicity, we consider point-like creases, although the model can be generalized to extended creases where a crease is a localized regions with a given natural curvature, as shown in reference Jules et al. (2019). The variation of functional (3) yields a set of $n$ ordinary differential equations given by (see Appendix A)

[TABLE]

where $\{c_{i}\}^{n}_{i=1}$ is the set of $n$ integration constants. The above equation describes the equilibrium shapes of the Euler’s Elastica. In the present work we assume that all the panels have the same constant, namely, $c_{i}=c$ for $i=1,\dots,n$. By comparing with the linear model of $f$-cones, one notices that $-c$ is proportional to the hoop stress $\sigma_{\varphi\varphi}$ in the limit of small deflections, thus, the sign of $c$ dictates whether the structure is in azimuthal compression ($c>0$) or tension ($c<0$) Lechenault and Adda-Bedia (2015); Cerda and Mahadevan (2005). The hypothesis of equal constants holds provided that there are no external forces acting on the creases introducing additional stresses in different panels, so that $\sigma_{\varphi\varphi}$ is continuous across the panels. From varying Eq. (3) one also obtains boundary terms that combine with terms coming from the variation of the energy stored in the creases. These boundary terms give the natural boundary conditions to solve equation (4) for each panel. In absence of external forces, these terms must satisfy

[TABLE]

where $\mathbf{f}_{i}=\kappa^{\prime}_{i}\mathbf{n}_{i}-\left(\kappa^{2}_{i}/2+c_{i}\right)\mathbf{t}_{i}$, which can be interpreted as a normalized force per unit-length along a ray of fixed $r$ Guven and Müller (2008).

At this point, it is convenient to introduce the vector $\mathbf{J}\equiv-\mathbf{u}\times\mathbf{f}+\kappa\mathbf{u}$ which is a conserved quantity associated with the rotational invariance of the system and can be interpreted as a torque Guven and Müller (2008). It can be shown that the quantity $J^{2}-c^{2}$ corresponds to the first integral of equation (4). One can use the vector $\mathbf{J}$ to obtain the equilibrium shape of the $f$-cone by first setting it parallel to the $z$-axis and projecting it onto the frame $(\mathbf{u},\mathbf{n})$, obtaining $\mathbf{J}\cdot\mathbf{u}=J\cos\theta=\kappa$ and $\mathbf{J}\cdot\mathbf{n}=J\varphi^{\prime}\sin^{2}\theta=\kappa^{2}/2+c$.

II.1 Infinitely stiff creases

In this section, we solve the case of infinitely stiff creases, so that $\delta g_{i}=0$. This means that the set of folding angles $\{\psi_{i}\}^{n}_{i=1}$ is an input of the problem and that the final solutions are parameterized by these angles. In Appendix B, we show that the boundary terms (5), together with the condition $\delta\psi_{i}=0$, imply the following:

[TABLE]

which means that the curvature is continuous through the crease. Also, the transversal force $\mathbf{f}$ is continuous, which can be written as

[TABLE]

The continuity conditions (6,7) imply that $\mathbf{J}^{+}_{i}=\leavevmode\nobreak\ \mathbf{J}^{-}_{i}$, thus, the entire structure is characterized by a single vector $\mathbf{J}$.

Solving the system (4) with the assumption of equal constants, $c_{i}=c$, requires $2n+1$ boundary conditions. Combining equations (6) and (7), one can show that

[TABLE]

thus, yielding $2n$ boundary conditions. Adding the closure condition (2) makes the problem well-posed.

Integration of Eq. (4) for each panel gives two possible solutions:

[TABLE]

where $\text{cn}(\cdot|\cdot)$ and $\text{dn}(\cdot|\cdot)$ are the cosine and delta amplitude Jacobian elliptic functions Abramowitz and Stegun (1965), with parameters $m^{\mbox{\tiny R}}_{i}$ and $m^{\mbox{\tiny S}}_{i}$ given by

[TABLE]

Here, $\kappa_{0i}$, $S_{0i}$ and $c$ are $2n+1$ unknown parameters that must be fixed such that the boundary conditions and the closure conditions are satisfied. Without loss of generality, one can define the cosine and delta amplitude functions such that $0<m^{\mbox{\tiny S}}_{i}<1$ and $0<m^{\mbox{\tiny R}}_{i}<1$ Abramowitz and Stegun (1965). Therefore, Eq. (11) shows that there exists two families of solutions. We identify as the rest configuration (thus the superscript R) the solutions for which $J^{2}>c^{2}$ and the snapped configuration (superscript S) with $J^{2}<c^{2}$. This choice is in agreement with the fact that we find a posteriori that the bending energy of the rest configuration is the lower one. As the quantity $J^{2}-c^{2}$ is defined for the entire structure, all the panels will be in either a rest state, or a snapped state, therefore, there is no mixture of states (provided that there is a single constant $c$). Notice that the present nonlinear approach allows us to justify the existence of two families of solutions for a given set of folding angles, regardless of whether the cone is symmetrical or not, a property that is difficult to prove in the linear model.

From this point forward, we consider a simplified situation of symmetrical $f$-cones made of all-mountain creases with same folding angle $\psi$ (we shall also omit the subscripts $i$ labelling the panels). Thus, it is sufficient to solve Eq. (4) in a single panel where $-\alpha/2\leq s\leq\alpha/2$ and $\alpha=2\pi/n$. For both rest and snapped states, we proceed to find numerically the parameters $\kappa_{0}$, $S_{0}$ and $c$ that satisfy the closure condition (2) and the boundary conditions (8,9). The closure condition (2) reads now $\Delta\varphi=\pm\alpha$ or [math]. For each state, the quantity $\Delta\varphi$ is an integral that can be computed analytically (see Appendix C). Then, for a given $S_{0}$, the closure condition defines a curve in the parameter space $\{\kappa_{0},c\}$ for each state. As the $f$-cone is symmetrical, the solutions given by Eq. (10) should be even functions with respect to $s=0$. Therefore, one has $S_{0}=0$ for the rest state and, for the snapped state there are two possibilies, either $S_{0}=0$ or $S_{0}=K(m_{s})$ ($K(\cdot)$ is the complete elliptic integral of the first kind and corresponds to the half period of the function $\text{dn}(\cdot|m)$) Abramowitz and Stegun (1965). We only found solutions for $S_{0}=K(m_{s})$ that satisfy the closure condition. Eq. (8) is automatically satisfied by the symmetry of the solutions. Then, Eq. (9) defines a second curve in the parameter space $\{\kappa_{0},c\}$ (one for each state). Hence, the solution for a given folding angle $\psi$ corresponds to the intersection between these two curves, so that the values of $\kappa_{0}(\psi)$ and $c(\psi)$ are obtained. This procedure is done for all $\psi\in[0,\pi]$ ($\psi>\pi$ would correspond to equivalent states but vertically flipped).

Hereinafter, we denote as $n^{\mbox{\tiny R}}$ (resp., $n^{\mbox{\tiny S}}$) the all-mountain $f$-cones with $n$-creases in a rest (resp., snapped) state. We show here the solutions for the most relevant cases: a semi-infinite crease ($1^{\mbox{\tiny R},\mbox{\tiny S}}$), a single infinite crease ($2^{\mbox{\tiny R},\mbox{\tiny S}}$) and two perpendicular mountain creases ($4^{\mbox{\tiny R},\mbox{\tiny S}}$). Fig. 3 shows the resulting $c(\psi)$ for these different cases. We noticed that except for the trivial case $2^{\mbox{\tiny R}}$, one has $c\nrightarrow 0$ as $\psi\rightarrow\pi$, suggesting the existence of a residual hoop stress as one approaches the flat state. Indeed, the residual hoop stress at $\psi=\pi$ coincides with the values of the hoop stress in the linear model Lechenault and Adda-Bedia (2015). We attribute this residual stress to a critical load needed to observe buckling of the facets, as it happens in Euler-Bernoulli beam buckling. It has been found that the Lagrange multiplier associated with the hoop stress has a fixed value for each state. The Lagrange multiplier $c$ depends on the folding angle, except for $1^{\mbox{\tiny S}}$ where $c(\psi)$ is nearly constant. In the configuration $4^{\mbox{\tiny R}}$, $c$ changes sign as $\psi\rightarrow 0$, suggesting at first sight a modification of the hoop stress from compressive to tensile as the folding angle gets sharper. However, this could be misleading because the interpretation of $c$ as a hoop stress is valid only for small deformations.

The equilibrium shapes of the $f$-cones are plotted using the coordinates $\theta(s)$, $\varphi(s)$. Thus, they are universal in the sense that they depend only on the folding angle $\psi$ and not the material properties. Nevertheless, changing the bending modulus or the dimensions of the cone will modify the stresses and torques supported by the system as well as its elastic response. Fig. 4 shows the deviations of the structure from the flat state. To quantify such deviations, we use the angle $\beta=\pi-\cos^{-1}(\mathbf{u}(0)\cdot\mathbf{u}(\pi))$ for a $1^{\mbox{\tiny R},\mbox{\tiny S}}$ $f$-cone and the polar angle $\theta_{c}=\cos^{-1}(\mathbf{z}\cdot\mathbf{u}(\alpha/2))$ for $n^{\mbox{\tiny R},\mbox{\tiny S}}$ $f$-cones ($n>1$). The choice of $\psi$, instead of $\theta_{c}$ (or $\beta$ for a $1^{\mbox{\tiny R},\mbox{\tiny S}}$ $f$-cone) as a control parameter prevents finding unphysical self-intersecting solutions. In Fig. 4(c), we compare $\theta_{c}(\psi)$ for $2^{\mbox{\tiny S}}$ configuration with the prediction of the linear model $\theta_{c}=\pi/2+0.4386(\psi-\pi)$ Lechenault and Adda-Bedia (2015); Walker and Seffen (2018). The polar angle of the crease deviates slightly from the linear prediction for sharper folding angles. This deviation fits better the experimental results shown in reference Walker and Seffen (2018), although there is still an important discrepancy from the theory because of the limits of the inextensibility hypothesis, which we explore in more detail in section IV.

II.2 Crease mechanics

The mechanical energy of a single point-like crease can be written as a function of invariants built from the three unit vectors that define the crease geometry: the crease vector and two vectors tangent to each facet Brunck et al. (2016). Such invariants are $\mathbf{t}^{-}_{i}\cdot\mathbf{t}^{+}_{i}$ and $\left(\mathbf{t}^{-}_{i}\times\mathbf{t}^{+}_{i}\right)\cdot\mathbf{u}^{+}_{i}$. At leading order, the elastic energy of the $i$-th crease takes the form Brunck et al. (2016)

[TABLE]

where $L=R-r_{0}$, $\sigma_{i}$ and $\tau_{i}$ are material constants associated to the crease. The crease energy can be rewritten in terms of the folding angle $\psi_{i}$, defined as the oriented angle $(\widehat{-\mathbf{t}^{-}_{i},\mathbf{t}^{+}_{i}})$ (see Fig. 2). Introducing the constants $k_{i}$ and $\psi^{0}_{i}$, such that

[TABLE]

allows us to rewrite (12) as $g_{i}(\psi_{i})=-L\,k_{i}\,\cos\left(\psi_{i}-\psi^{0}_{i}\right)$. Thus, the crease energy $g_{i}=g_{i}(\psi_{i})$ is a function of the folding angle which is now an unknown variable. If $\psi_{i}\approx\psi^{0}_{i}$, the crease energy approximates to the energy of an elastic hinge $g_{i}\approx\leavevmode\nobreak\ L\,k_{i}\,(\psi_{i}-\psi^{0}_{i})^{2}/2+E_{0}$, where $E_{0}$ is a constant, $k_{i}$ is the crease stiffness, and $\psi^{0}_{i}$ is the rest angle of the crease. The conical geometry implies that the folding angle is constant along the crease. One can show that taking into account the terms coming from the variation of the crease energy, we obtain an additional boundary condition (see Appendix B)

[TABLE]

Eq. (14) states that the value of the curvature at the crease is given by the moment imposed by its mechanical response. For a crease energy given by Eq. (12), one has $\kappa^{+}_{i}=\bar{k}_{i}\,\sin\left(\psi_{i}-\psi^{0}_{i}\right)$, where $\bar{k}_{i}=Lk_{i}/a$ is the normalized crease stiffness.

We study again configurations with equally spaced mountain creases that have the same mechanical properties, such that, $\bar{k}_{i}=\bar{k}$ and $\psi_{i}^{0}=\psi_{0}$ . The control parameters are then the crease stiffness $\bar{k}$ and the rest angle $\psi_{0}$. The boundary conditions are those of the infinitely stiff creases case supplemented by Eq. (14). Therefore, one can use the solutions found for the infinitely stiff creases and search the value of $\psi$ such that Eq. (14) is satisfied.

Figures 5(a-c) show typical shapes for the cases $1^{\mbox{\tiny R},\mbox{\tiny S}},2^{\mbox{\tiny R},\mbox{\tiny S}},4^{\mbox{\tiny R},\mbox{\tiny S}}$. Fig. 5(d) shows the final folding angle $\psi$ as function of $\bar{k}$ for a fixed rest angle $\psi_{0}$ of the crease. Notice that $\psi\rightarrow\pi$ as $\bar{k}\rightarrow 0$ and that $\psi\rightarrow\psi_{0}$ as $\bar{k}\rightarrow\infty$. For a given $\bar{k}$, the snapped state always displays a larger $\psi$ than the rest state. This observation is explained by the fact that the hoop stress of the snapped state is always larger than its respective rest state.

The total energy of the structure can be computed by summing the bending energy of the facets and the energy of the creases. Fig. 6(a) shows an example of the energy landscape as function of the polar angle $\theta_{c}$ of the crease for an all-mountain $f$-cone with 4 creases. While the bending energy of the facets has an asymmetric parabolic shape, the energy of the creases has a double-well potential shape whose minima are at the same energy level. The resulting shape of the total energy is an asymmetric double-well potential where the two minima corresponds to the rest and the snap states. Figure 6 (b) shows the derivative of the energy with respect to $\theta_{c}$ which corresponds to the moment $M(\theta_{c})$ applied on the creases. Notice that the energy has a cusp at $\theta_{c}=\pi/2$ which leads to a discontinuity in the mechanical response of the structure. This type of snap-through transition has been found in similar systems such as the waterbomb origami with rigid facets Hanna et al. (2014). However, the snap-through mechanisms of the $f$-cone and the waterbomb are different. The former is mediated by the asymmetry in the bending energy between the rest and snapped states and the latter by the asymmetrical kinematical conditions imposed by the rigid facets. For the waterbomb, the only relevant energy is the mechanical energy stored in the creases which, in addition to the kinematical constraints, accounts for an asymmetric double-well like potential.

III Continuous elastic model of origami structures

Commonly, the mechanical response of origami-based metamaterials can not be reduced to elastic hinges connected to rigid or isometric panels Silverberg et al. (2015). The energy landscape of deformations depends generally on both bending and stretching energies of the panels as well as on the inherent spatially extended nature of the creases Jules et al. (2019). Therefore, one needs to supplement the present analysis with a more accurate description that takes into account these different contributions. When applied to the $f$-cone, such a description should be validated by the bounds given by the analytical model.

In the following, we propose a mechanical model based on a continuous description of origami structures that is suitable for numerical implementation and test it on the $f$-cone by performing Finite Element analysis (FEA). The simulations were designed in the commercial FEA software COMSOL Multiphysics 5.4. Within this software, the Structural Mechanics Module is equipped with quadratic shell elements, which have been used to our purpose. All the simulations were carried out with a linear elastic Hookean material model and geometric non-linear kinematic relations have been included. The plate Young’s modulus is $E=3.5$ GPa, the Poisson’s ratio is $\nu=0.39$, and the plate thickness is $h=300\,\mu$m. We searched for solutions with the default stationary solver, where the non-linear Newton method has been implemented. Mesh refinement studies were undertaken to ensure convergence of the results.

III.1 Temperature-induced hingelike creases

In our numerical model, we develop a method to create creases that are able to reproduce the hingelike mechanical response of commonly folded thin sheets. Inspired by the experimental results of Ref. Jules et al. (2019), we model creases as narrow slices of the plate undergoing thermal expansion due to a temperature gradient through the thickness of the plate, as it is schematically shown in Fig. 7.

In order to test the mechanical response of these temperature driven creases, we first perform a numerical test of a single fold in a hingelike geometry consisting of two facets and the crease in the middle. We consider a rectangular plate of length $L$, width $W$ and thickness $h$. In the middle of the plate, we take a transversal narrow slice of width $b\ll W$, across the width of the plate, dividing the plate into three parts. The narrow slice corresponds to the crease while the other two parts correspond to the facets. A linear temperature gradient $\Delta T$ is applied through the thickness of the plate. We let the creased region of the plate to undergo thermal expansion by defining an inhomogeneous coefficient of thermal expansion $\alpha_{T}(s)$ which is constant for $|s|<b$ and zero elsewhere, where $s$ is the arc-length perpendicular to the creases. To simulate a more realistic crease, we add a rigid connector made of an infinitely stiff material of length $W$, thus preventing bending deformation in the direction of the crease line. Moreover, tuning the mechanical response of the crease is done by varying the thickness $h_{c}$ and the Young’s modulus $E_{c}$ of the crease slice. The parameters $h_{c}$, $E_{c}$, $\Delta T$ and $b$ will define the crease mechanical response.

Taking advantage of the two-plane symmetry of the single fold geometry, we solve only for one quarter of the plate and then obtain the entire equilibrium shape by reflections. We perform two different studies: a heating up test and a mechanical response test. The first study consists in heating the plate from the bottom while the ends of the two plates parallel to the crease are constrained to move in the $xy$-plane. When heating, the crease bends towards the sense of lower temperature (bottom), while the facets remain practically flat. The resulting angle between the two facets corresponds to the rest angle $\psi_{0}$ of the crease. The second study consists in a test of the mechanical response of the crease. We add four additional rigid connectors to the sides of the facets that are perpendicular to the fold line. Two moments of opposite signs are applied respectively to each pair of rigid connectors attached to the facets so that the hingelike system can open or close. Through this mechanical test, we were able to verify the hingelike behavior of the crease.

In the following, we employ the temperature-induced creases in our numerical model of $f$-cones.

IV Numerical analysis of foldable cones

We begin with a circular planar disc of external radius $R=100$ mm and a central hole of radius $r_{0}=1$ mm. Then, $n$ radial narrow slices of constant width $b$, that correspond to the creases of the $f$-cone, are created. Rigid connectors along the creases are added so as to prevent bending along the longitudinal direction of the creases. In order to take advantage of the symmetries of the system, depending on the number of creases, only the fundamental unit cells are numerically solved and then the complete structure is reconstructed using reflections through the symmetry planes. Because each plane of symmetry cannot coincide with a rigid connector, the $1^{\mbox{\tiny R},\mbox{\tiny S}}$ case must be solved entirely. For $2^{\mbox{\tiny R},\mbox{\tiny S}}$, only a half of the disk is numerically solved so that the axis of symmetry is perpendicular to both creases. For $4^{\mbox{\tiny R},\mbox{\tiny S}}$, a quarter of the disk is solved so that a single crease is at $45\textordmasculine$ from one plane of symmetry.

IV.1 Indentation tests

Hereinafter, we focus only in the all-mountain $f$-cone with 4 creases, anticipating that our general conclusions also apply to more complex configurations. In order to test our analytical predictions, we study the snapping of the system through an indentation process from the rest state to the snapped state, which is carried out in two steps. Initially, we turn on the temperature to take the $f$-cone to its rest state with a given folding angle $\psi$. The ends of the creases are constrained to move in the $xy$-plane, so that the points at the central rim rise up when the temperature is activated. For all our simulations, we fix the temperature gradient in the crease such that $\alpha_{T}\Delta T=0.2$. We also choose $b=1$ mm, which is within the expected order of magnitude with respect to $h$ according to crease formation measurements in thin sheets Benusiglio et al. (2012).

In a second step, the central rim is vertically lowered quasistatically while constraining the ends of the creases to move in the $xy$-plane. Each crease is constrained to rotate only in the plane defined by its initial direction and the $z$-axis (as shown in Fig. 4(a)). Throughout indentation, the vertical displacement of the central rim is specified and a reaction moment $M(\theta_{c})$ at the creases is computed. Therefore, the mechanical response of the $f$-cone to the indentation process consists in the determination of the curve $M(\theta_{c})$ (see movies available as supplemental material of such indentation tests in Sup ).

The folding angle $\psi$ of the creases is also tracked during indentation. To measure $\psi$ from the numerical results, we extract the tangent vector field of concentric curves initially defined in the flat configuration. We evaluate this vector field at each side of the crease and measure the resulting angle between them. The local folding angle is found to be not exactly constant along the crease but is a function of the radial coordinate. For this reason, a representative measurement of $\psi$ is chosen to be the average between two local folding angles measured at radial distances $r_{0}+(R-r_{0})/3$ and $r_{0}+2(R-r_{0})/3)$. In each indentation test, we extract a curve $\theta_{c}$ as function of $\psi$ which we call the indentation path.

The resulting curves $M(\theta_{c})$ and $\theta_{c}(\psi)$ will be discussed in the following.

IV.2 Numerical results

Fig. 8 shows the variation of $\theta_{c}(\psi)$ obtained numerically during the indentation test for different crease thicknesses and Young moduli. These results are compared with the polar angle $\theta_{c}(\psi)$ given by analytical $f$-cone calculations shown in Fig. 4(d). This parametric study allows us to highlight to what extent the softness of the crease affects the indentation path. Our results show that our continuous model of the crease is more sensitive to variations of $h_{c}$ than those of $E_{c}$. We notice that the indentation paths do not generally follow the analytical solutions given by the isometric constraint, meaning that the intermediate shapes throughout indentation are not perfect developable cones. If the crease is too stiff, as in the case of Fig. 8(a), the indentation path follows a nearly vertical line (i.e. approximately constant folding angle path) connecting the two stable points. However, when the crease is made softer (Fig. 8(b)), either by reducing its thickness or its Young’s modulus, the indentation paths approach the one predicted by the isometric constraint.

While the rest states are generally well predicted by the isometric constraint, the snapped states depart from the analytical predictions when $h_{c}$ is decreased. This result could be explained in terms of the crease stiffness. For a stiff crease, the folding angle is roughly constant along the crease enforcing the shape to be closer to a perfectly developable cone. On the other hand, a $4^{\mbox{\tiny S}}$ cone is characterized by an azimuthal tension which favors stretching deformations of the panels and thus tensile traction on the crease. This mechanism could induce a varying folding angle along the crease and yield a structure that can depart from a perfect developable cone, especially for softer creases. To verify this analysis, we plot in Fig. 9 one quadrant containing the lines of smallest principal curvature for stiff and soft crease cases. In a perfect developable cone, these lines coincide with the generators of a surface (i.e. lines of zero curvature), however, we observe that the lines curve significantly close to the vertex, where the energy is concentrated. This effect is more pronounced for the soft crease case than the stiff one.

Fig. 10(a) shows typical curves for the moment as function of the polar angle $\theta_{c}$ during indentation for both stiff and soft crease. Notice that the unstable branch of the stiffer crease is higher than that of the soft one, which in a real experiment leads to more energy being released during a snapping process. This observation can be attributed to a larger stretching energy barrier that is required to be overcome, which is evident from the energy plots shown in Fig. 10(b). To compare with analytical predictions, one should focus on the Inset of Fig. 10(b) which shows the evolution of the bending energy of the facets throughout indentation and does not take into account the bending energy at the creases. The predicted bending energy exhibits an asymmetric parabolic shape, where the minimum corresponds to a flat solution. It is obvious that the bending energy of a soft crease follows closer the prediction than that of the stiff crease, which has a pronounced convex shape in the middle.

V Conclusion

Foldable cones are the simplest example of a single-vertex origami whose facets can bend. In the present work, we developed a theoretical model which allows us to obtain the shape of $f$-cones for any deflection. The model shows that the bistable behavior of these structures is robust, regardless the specific properties of the creases. In particular, for symmetrical all-mountain $f$-cones we obtained the polar angle at the crease as function of folding angle for both rest and snapped states.

However, in more realistic situations, the geometry and mechanical response of a $f$-cone are characterized by a competition between the elasticity of the facets (both their bending and stretching behavior) and the stiffness of the creases. To this purpose, we have developed a continuous numerical model accounting for the both the elasticity of the creases and facets. Applying to the particular case of two perpendicular mountain creases we numerically studied the role of crease stiffness and verified the snap-thorough behavior through a series of indentation tests. We studied the indentation paths in the $\theta_{c}(\psi)$ diagram and showed that the structures do not follow the shape of a perfect cone throughout the indentation. For stiff creases, the path followed is that of an approximately constant folding angle while the two stable states lie closely to the theoretical prediction. When the crease is made soft, the indentation paths follow closely the branches given by the isometrical constraints, however, it is noted that while the shape of the rest state is close to the theoretical prediction, the snapped one deviates further from it. From an energetic viewpoint, not only do stiffer creases lead to indentation paths with higher stretching energy barriers, but they also enforce the preferred angle more strongly. Hence, it can be concluded that a $f$-cone made with stiffer creases requires more stretching when passing through $\theta_{c}\sim\pi/2$. On the other hand, while softer creases induce large deviations of the preferred angles, they allow for low stretching during the inversion, which explains why they follow the bounds set by the analytical calculations more closely.

The present study validates our numerical model of temperature-induced hinge-like creases which can be applied to origami structures with more complex extended networks. In this case, a temperature-induced folding of the network would work as a phase-field model where the sharp piecewise energy landscape is replaced by a smooth curve. The choice of temperature field as a trigger for crease formation is arbitrary as any other diffusion field, such as concentration Klein et al. (2007) or swelling Kim et al. (2012), would play a similar role. The main sought mechanism is to build up a reference configuration with a noneuclidean reference metric due to the presence of an initially imprinted crease network Wu et al. (2013). This approach is advantageous since one does not need to track sharp boundaries where the deformation fields are discontinuous. It renders numerical implementation tractable and less time-consuming, two important aspects when implementing the mechanical behavior of complex origami or crumpled structures.

Acknowledgements.

I.A.-S. acknowledges the financial support of CONICYT DOCTORADO BECAS CHILE 2016-72170417. M.A.D. thanks the Velux Foundations for the support under the Villum Experiment program (project number 00023059).

Appendix A Derivation of the Euler Elastica Equation

The derivation of Euler’s Elastica from the energy functional given by Eq. (3) can be found in Guven and Müller (2008). Is it instructive to repeat the calculations here for the sake of completeness. The variation of the functional (3) gives

[TABLE]

where we have used the identities: $\mathbf{u}\cdot\mathbf{t}\times\delta\mathbf{t}^{\prime}=\mathbf{n}\cdot\delta\mathbf{t}^{\prime}$ and $\mathbf{u}\cdot\delta\mathbf{t}\times\mathbf{t}^{\prime}=-\kappa\mathbf{t}\cdot\delta\mathbf{t}$. The term $\delta g_{i}$ contributes to boundary terms only, and will be treated below. Taking $\mathbf{u}_{i}$ and $\mathbf{t}_{i}$ as independent variables, the terms proportional to $\delta\mathbf{u}_{i}$ yield

[TABLE]

while the terms proportional to $\delta\mathbf{t}_{i}$ give

[TABLE]

Notice that

[TABLE]

Differentiating Eq. (17) with respect to $s$ and using (18), it follows that $\Lambda^{\prime}_{i}=5\kappa_{i}\kappa^{\prime}_{i}$. Integrating once, we obtain $\Lambda_{i}=5\kappa_{i}^{2}/2+c_{i}$, where $c_{i}$ is an integration constant. Then, Eq. (17) can be written as follows

[TABLE]

Differentiating Eq. (20) with respect to $s$ and projecting onto $\mathbf{u}_{i}$ gives

[TABLE]

Projecting equation (16) onto $\mathbf{u}_{i}$ and equating with Eq. (21), one gets $\lambda_{i}=\kappa_{i}^{2}/2-c_{i}$. Then, Eq. (16) now reads

[TABLE]

On the other hand, one can diefferentiate once Eq. (20) and obtain

[TABLE]

Using equations (22,23), one obtains the Euler’s Elastica equations given by Eq. (4). In the following, the boundary terms will be treated.

Appendix B Boundary Conditions

It is useful to compute the variation of the functional (3) in terms of virtual rotations of the frame specified by the Euler-like angles. To this purpose, we first introduce the vectors $\mathbf{e}_{\varphi}=-\sin{\varphi}\,\mathbf{x}+\cos{\varphi}\,\mathbf{y}$ and $\mathbf{n}_{\varphi}\equiv\mathbf{u}\times\mathbf{e}_{\varphi}$ which span the plane containing the vectors $\mathbf{t}$ and $\mathbf{n}$ (for simplicity, we omit subscripts here). If $\phi(s)$ is the angle between $\mathbf{t}$ and $\mathbf{e}_{\varphi}$, then,

[TABLE]

Defining $\mathbf{e}_{\rho}=\cos{\varphi}\,\mathbf{x}+\sin{\varphi}\,\mathbf{y}$, one can write

[TABLE]

Then, one can show the following relations

[TABLE]

Notice that $\mathbf{u}\cdot\delta\mathbf{u}=0$ and $\mathbf{t}\cdot\delta\mathbf{t}=0$ as expected. The following relations are useful

[TABLE]

Now, we put the subscripts back and write some useful relations. First, notice that in the plane perpendicular to a crease, the frame $\{\mathbf{t},\mathbf{n}\}$ rotates by an angle $\psi_{i}-\pi$, which can be expressed as follows

[TABLE]

Using the relations given by Eq. (28) and Eq. (29) one obtains the following relations

[TABLE]

where $\phi^{\pm}_{i}$ is the angle between $\mathbf{t}^{\pm}_{i}$ and $\mathbf{e}_{\theta}$. Notice that $\psi_{i}=\pi+\phi^{+}_{i}-\phi^{-}_{i}$, then, $\delta\psi_{i}=\delta\phi^{+}_{i}-\delta\phi^{-}_{i}$. At this stage, we distinguish two cases: infinitely stiff creases and finite crease stiffness.

B.1 Infinitely stiff crease

By taking $\delta g_{i}=0$ and using the periodic convention in Eq. (5), one can write

[TABLE]

where we have used $\mathbf{u}^{-}_{i}=\leavevmode\nobreak\ \mathbf{u}^{+}_{i}$. Using Eqs. (30), the condition of infinitely stiff crease $\delta\psi_{i}=0$ is equivalent to imposing $\mathbf{n}^{+}_{i}\cdot\delta\mathbf{t}^{-}_{i}=\mathbf{n}^{-}_{i}\cdot\delta\mathbf{t}^{+}_{i}$. Thus, imposing $\delta\psi_{i}=0$ and letting $\delta\mathbf{u}^{+}_{i}$ undergo independent virtual rotation imply the boundary conditions (6) and (7). By projecting Eq. (7) onto $\mathbf{n}^{+}_{i}$ and $\mathbf{t}^{+}_{i}$, one obtains

[TABLE]

where we have assumed $c_{i}=c$. Manipulating Eq. (32) one obtains equations (8) and (9).

B.2 Finite crease stiffness

The variation of the crease energy given by Eq. (12) reads

[TABLE]

The last term in the right-hand side is zero because $(\mathbf{t}^{-}_{i}\times\mathbf{t}^{+}_{i})|_{s_{i}}$ is parallel to $\mathbf{u}_{i}(s_{i})$. Using the cyclic properties of the triple product we can write

[TABLE]

Using the identities (30), Eq. (34) can be rewritten as follows

[TABLE]

Notice that the above equation has the form $\delta g_{i}=(dg_{i}/d\psi_{i})\delta\psi_{i}$.

Using the definition for the constants introduced in Eq. (13) and recalling that $\theta^{+}_{i}=\theta^{-}_{i}$ and $\varphi^{+}_{i}=\varphi^{-}_{i}$, we can rewrite Eq. (5) as follows

[TABLE]

The infinitesimal variations of the frame vectors can be translated to virtual rotations in terms of the Euler angles $\delta\theta^{+}_{i}$, $\delta\varphi^{+}_{i}$, $\delta\phi^{+}_{i}$ and $\delta\phi^{-}_{i}$. Assuming that all these virtual rotations are independent, one obtains conditions (6), (7) and (14) with $g_{i}(\psi_{i})$ given by Eq. (12).

Appendix C Closure condition

In a symmetrical $f$-cone with $n$ creases, the azimuthal angle spanned by a single panel is given by the integral

[TABLE]

where $\alpha=2\pi/n$. As the deformed state will also be symmetrical, the closure condition can be written as $\Delta\varphi=\pm\alpha$ or $\Delta\varphi=0$ according to Eq. (2). These integrals can be computed analytically for each state:

[TABLE]

and

[TABLE]

where the labels $\mbox{R},\mbox{S}$ stand, respectively, for the rest and snapped states. Also, $\Pi(\cdot,\cdot|m)$ is the elliptic integral of third kind and $\text{am}(\cdot|m)$ is the Jacobi amplitude with modulus $m$ Abramowitz and Stegun (1965).