title: Cutting holes in bistable folds
authors: T. Yu I. Andrade-Silva M. A. Dias J. A. Hanna
text: 
  	
  	
  	tiany@princeton.edu
  	Department of Civil and Environmental Engineering, 
  Princeton University, Princeton, NJ 08544, USA
  	
  	
  	ignacio.andrade-silva@univ-amu.fr
  	Aix-Marseille University, CNRS, IUSTI, Marseille 13013, France
  		
  	
  	marcelo.dias@ed.ac.uk
  	Institute for Infrastructure & Environment, School of Engineering, The University of Edinburgh, Edinburgh EH9 3FG, UK 
  
  	Department of Mechanical and Production Engineering, Aarhus University, 8000 Aarhus C
  	
  	
  	jhanna@unr.edu
  	Department of Mechanical Engineering, University of Nevada, 1664 N. Virginia St. (0312),  Reno, NV 89557-0312, USA
  	
  
  A folded disk is bistable, as it can be popped through to an inverted state with elastic energy localized in a small, highly-deformed region on the fold. Cutting out this singularity relaxes the surrounding material and leads to a loss of bistability when the hole dimensions reach a critical size.  These dimensions are strongly anisotropic and feature a surprising re-entrant behavior, such that removal of additional material can re-stabilize the inverted state.  A model of the surface as a wide annular developable strip is found to capture the qualitative observations in experiments and simulations.  These phenomena are consequential to the mechanics and design of crumpled elastic sheets, developable surfaces, origami and kirigami, and other deployable and compliant structures.
  
  Cutting holes in bistable folds
      J. A. Hanna
      December 30, 2023
  ===============================
  
  
  
  
  The role of elastic singularities in the deformation of thin sheets and shells is still poorly understood, despite a quarter of a century of intense investigation into their geometry and energetics <cit.>.
  Over the years, several perspectives have emerged, viewing these localized high-energy regions as a manifestation of spontaneous condensation of both curvature and stretching  
  <cit.>,
  sources of rigidity <cit.>, 
  or an organizing framework for random crumpling <cit.>,
  regular patterns <cit.>,
  or dynamics <cit.>.
  Our interest in the current study was driven by Witten's observation <cit.> that the excision of such regions of focused elastic energy leads to significant relaxation of neighboring regions of material, and the indication that such surgery should also modify the rigidity and stability landscape of any surrounding structure.
  The proximity of the edge of regression, or other virtual singularities living outside nominally inextensible surfaces, has been qualitatively linked to the structural stiffness response <cit.>.  
  Another thread in this work is the question of multistability of systems of creases and facets with competing flexibilities, including origamic analogs of elastic singularities <cit.>. 
  
  Our model system is a single fold in an elastic disk, and the singular structure formed by popping it through with a thumb.  
  This is perhaps the simplest bistable “foldable cone” examined in <cit.>.  
  In <cit.> it was noted that a small hole reduced the energy barrier to pop through a fold to its inverted state, but this was not pursued to its logical conclusion, the complete elimination of the barrier with a sufficiently large hole.
  In this note, we employ theory, numerics, and experiment to capture the complex behavior of a fold after removal of its singularity and a variable zone of surrounding material.
  We find that cutting a hole of sufficient size around the singularity leads to a loss of bistability through a fold bifurcation that destabilizes the inverted state.  There is significant anisotropy in the critical hole dimensions, such that a narrow slit aligned along the crease can be as large as the disk without destroying bistability.  We also observe a curious re-entrant behavior of the stability diagram for small elliptical holes aligned perpendicular to the crease, which can in some parameter ranges more effectively eliminate bistable behavior than a larger circular hole. 
  We demonstrate the surprising applicability of a developable ribbon model to this class of wide, topologically annular shapes.
  
  Experiments were performed on disks of radius 1 ≤ R ≤ 10 cm, thickness t= 0.005 in (0.127 mm), 0.003 in (0.076 mm), or 0.002 in (0.051 mm), and prescribed central elliptical hole geometries (semiaxes a and b perpendicular and parallel to the eventual crease), obtained by cutting (Cameo 3, Silhouette America, Lindon, UT) polyester shim stock (Artus Corp., Englewood, NJ) and subsequently creasing along a diameter using a vise set to a prescribed position and held for two seconds.
  The structure was then flexed by inverting it once or twice and allowed to relax for five minutes before bistability tests were performed, with the sample hanging such that the crease was vertical to minimize gravitational effects on the bistability of the thinnest sheets.
   We refer to the simply creased state as the folded configuration 
    and the stable popped-through state as the inverted configuration. 
   When both stable states exist, they are separated by an energy barrier corresponding to another, unstable, equilibrium state.  This barrier was examined in <cit.> for a reduced model of the surface with a single hole size. 
  The folded state is characterized by a rest crease angle γ_0 over which we have little control; thickness, disk size, and hole size all contribute significantly to this value, which we report as a range spanning multiple hole sizes for a given sample set (measured from photographs of the disks).
  The inverted state is characterized by a final crease angle γ_f, which is observed to be a function of radial position in experiments and numerics, and an angle δ between a side of the crease and a line connecting the ends of the two sides.  Figure <ref>(a-b) illustrates examples of the two states and associated parameters for two different hole geometries.
  The inverted configuration reflects a competition between the unknown and uncontrollable stiffness of the crease and the bending resistance of the facets comprising the remainder of the disk, such that the overall disk size R is a relevant scale that we can understand using the related concepts of “origami length” <cit.> or “hinge index” <cit.>.
  We should expect an asymptotic approach to a linear scaling of critical bistable hole size with the disk radius as the latter grows and the crease becomes effectively rigid.  However, the interference of gravity also becomes more important with increasing disk radius, setting a practical limit for the experiments.
  
  
  
  We employ an annular wide-strip model that treats the punctured disk as a developable surface outside of the crease, an approach that is reasonable for static configurations of sufficiently thin elastic sheets <cit.>.
  In this model, we treat the crease as a generator (zero-curvature direction) with uniform final angle γ_f, an approximation we will revisit shortly below.
  Full details of the model, including boundary conditions and numerical implementation, can be found in Appendix <ref>; we sketch the important aspects here.
   The directrix r(s) forms the outer circumference of one half of the symmetric disk, parameterized by arc length 0 ≤ s ≤π R, and carries an orthonormal material (Darboux) frame of curve tangent T = r', where a prime denotes an s-derivative, surface normal N, and surface tangent normal B = T×N.
   The evolution of the frame is given by T'= κ_n N  - κ_g B, N'= -κ_n T + τ_g B, B'= κ_g T - τ_g N, where  κ_n, κ_g = -1/R, and τ_g are the normal curvature, geodesic curvature, and geodesic torsion.
    As shown in Figure <ref>(c), generators lie in the B+ηT direction, making a local angle β with the directrix; η = β = τ_g / κ_n.
   The shape is symmetric and given by the embedding
  
      X(s,v)=r(s)+v[B(s)+η(s)T(s)]   ,
  
  with 0 ≤ v ≤ V(s,η  ; R,a,b) the coordinate along the generator of implicitly treated length V √(1+η ^2).
  This surface has mean curvature H= κ_n (1+η^2)2[1+v[η'+κ_g(1+η^2)]] and  area element dA=[1+v(η'+κ_g (1+η^2))]   dsdv. 
  Defining a crease stiffness per unit length K_c and a facet bending rigidity D=E t^3 /[12 (1 -ν ^2)] incorporating the Young's modulus E and Poisson's ratio ν, the total elastic energy U of a creased punctured disk can be written as 
  an augmented Wunderlich functional <cit.>
  with contributions from both crease and facets,
  
  
  
      U/2D   = K_c/D(R-b)  ∫_γ_0^γ_fsin(γ̃_f-γ_0) d γ̃_f 
       
        + 1/2∫_0^π R∫_0^V(2H)^2  dA   ,  
         = K_c R/D(1-b/R) [1-cos (γ_f - γ_0)] + ∫_0^π R YW ds   ,
  
  
  with Y=κ^2_n (1+η^2)^2/2[η'+κ_g(1+η^2)] and W=ln[1+V (η'+κ_g (1+η^2))].
  The crease stiffness diverges as the origami length D/K_c approaches zero with the thickness <cit.>.
  The crease contribution to the energy (<ref>) only enters the problem through the boundary conditions, and does not appear in the Euler-Lagrange equations, given by <cit.>
  
      F'    = 0 , 
      M' + T×F   =0 , 
      ∂_κ_n (YW) -η M_1 -M_3    =0  , 
      ∂_η (YW) -(∂_η' (YW))'-κ_n M_1    =0  ,
  
  where forces and moments, normalized by D, have components in the moving frame specified by F=F_1 T+F_2 N+F_3 B and  M=M_1 T+M_2 N+M_3 B. 
  Equations (<ref>-<ref>) for one half of the symmetric structure, along with an Euler angle description of the moving frame, and boundary conditions imposed at the crease, are solved using the continuation package AUTO 07P <cit.>.
  
  To explore beyond the limitations of the developable model, as well as to allow independent control of rest crease angle, crease stiffness, and material thickness not possible in experiments, simulations were performed using the commercial finite element (FE) software COMSOL Multiphysics 5.4, employing quadratic shell elements, linear-elastic material, and creases introduced using a through-thickness thermal gradient <cit.>. 
  These simulations also avoid self-contact effects that are an issue with experiments at small hole sizes.  
  Full simulation details can be found in Appendix <ref>. 
  
  
  
  Results from the developable model (<ref>-<ref>) are shown in 
  Figure <ref>(a-c) for circular holes (a=b), rest crease angle γ_0=45^∘, and several values of crease stiffness K_cR/D.   
  Certain features are shared by all the solution curves in parameter space.  The stable inverted state and another unstable state are lost through a fold bifurcation at a critical hole size, which shows only moderate variation with crease stiffness.  At the end of the unstable lower branch, the curves terminate due to a failure of the embedding of the developable surface as the edge of regression approaches the sheet boundary.
  For stiffer creases, a greater proportion of the energy U is associated with facet bending energy U_b rather than the opening of the crease.
  At small hole sizes, the crease energy increases rapidly, and the facet energy does likewise for stiff creases but decreases for softer creases.
  Renderings of stable, near-fold-point, and unstable inverted states of the developable model with γ_0=45^∘ and K_cR/D=20 show the generators (black lines), edges of regression (red curves), and color maps of twice the squared mean curvature 2H^2 (with R set to unity) that appears in the facet energy.  Further details of these solutions can be found in Appendix <ref>.
  A small hole size leads to a shape similar to a generalized cone, which would have a single point “inside” the annulus as its edge of regression, rather than the cusped curve of the more general structure.  As the hole size increases, the generators near the crease first align with it, becoming “cylindrical” at a/R ≈ 0.116, echoing the qualitative features of the ridge relaxation observed in Witten's experiment <cit.>.  Past this point, the generators near the crease converge towards the outer perimeter of the annulus, so that a portion of the edge of regression lies “outside” it.
  On a developable strip, the bending moment is inversely proportional to the distance of a point on a generator to the edge of regression.  Thus, we would expect that if the developable constraint were relaxed to allow a non-uniform crease angle γ_f, 
  the inner part of the crease would open to a wider angle than the outer part for small holes and vice versa for large holes.  Indeed we observe this effect in experimental samples.  Results on this effect from a more quantitative analysis using FE are shown in 
  Figure <ref>(d-e) for the stable inverted configuration and the unstable state on the energy barrier, respectively, for several circular hole sizes up to a value close to the loss of bistability.  A representative set of crease parameters are used (the value R/20t can be considered an effective crease stiffness for the FE results, as discussed in Appendix <ref>). 
  In the inverted configuration, the crease angle variation along the normalized distance on the crease Y_C/R from the center of the hole shows the expected transition between inner and outer crease opening as the hole size increases.  The developable assumption of constant crease angle works better for large hole sizes.  In contrast, the unstable equilibrium is always more open towards the outer edge of the annulus.  Curiously, the crease angle near the hole in the unstable equilibrium is actually smaller than the rest angle.
  
  
  
  We next explore the anisotropic interaction between material removal and the elasticity of the structure by considering elliptical holes with semiaxes a perpendicular to and b parallel to the crease.
  Bistability boundaries for the developable model are shown in Figure <ref>(a) for several values of rest crease angle γ_0 and crease stiffness K_cR/D.
  Also shown are boundaries determined from experiments 
  on two thicknesses of material and from FE simulations of one thickness and set of crease parameters. A quantitative comparison across the three approaches is not possible. In experiments, the crease stiffness is unknown and the rest angle has a strong dependence on thickness and hole geometry; we report an average and standard deviation across a range of hole sizes in order to treat hole size as if it were an independent parameter.
  Thickness appears in the stiffness in the developable model; a previous observation indicates that the origami length D/K_c ≈ 200t for mylar (a polyester) sheets <cit.>, 
  
  although we note that the geometry of the current problem is significantly different than that of the cited work.  The crease becomes effectively rigid as the thickness vanishes.
  The FE simulations use an angle close to one of the developable solutions, with an effective stiffness falling within a range in which the dependence of the developable model on stiffness is very weak.
  However, the qualitative behavior of all the boundary curves is the same, and the quantitative spread is fairly small within a wide range of reasonable parameters.  Boundaries are shifted to larger hole sizes by a smaller crease angle (sharper fold) or by a stiffer crease; the experiments on thinner materials behave accordingly as having both sharper folds and stiffer creases. 
  An unexpected re-entrant behavior of the boundary curves is present in all three approaches: theory, experiments, and simulations.  This means that in some parameter ranges, removing more material actually leads to a reappearance of bistable behavior.  Fixing the length of the hole axis perpendicular to the crease, elliptical holes with long axis perpendicular or parallel to the crease can be monostable while less eccentric holes are bistable.  It is also apparent that long, slit-like holes along the crease do not eliminate bistability, while small perpendicular slits do.
  
  Solution curves, and renderings of stable inverted states near the bistability boundary, for the developable model with γ_0=45^∘ and K_cR/D=20, are shown in Figure <ref>(b-c).  
  Stable (upper) and unstable (lower) inverted states are seen to appear via an isola-center bifurcation.
  The edges of regression for elliptical hole inverted states can take more complicated multi-cusp forms not observed with circular holes.
  Further details of these solutions can be found in Appendix <ref>.
  
    
  
  In conclusion, we have examined the excision of high-energy material around an elastic singularity formed by inverting a simply folded thin disk. This process eliminates a source of rigidity, increasing the flexibility of the system. It reorients the low-curvature directions around the fold, influencing the opening angle distribution and eventually eliminating the inverted state, and thus the bistability, through a highly anisotropic mechanism re-entrant in the space of hole geometric parameters.
  These findings have consequences for the mechanical compliance and energetics of perforated thin sheets, and for the design of deployable structures, in which fatigue of a highly stressed vertex is undesirable, prompting the introduction of gaps. Beyond folded and cut structures, similar mechanics is expected in other compliant mechanisms featuring networks of hinges, facets, and springs.
  Bistability and critical hole dimensions are also influenced by the presence of multiple folds or the addition or removal of angular sectors of material to adjust the strength of the conical singularity, topics to be explored in detail in a future study <cit.>.
  
  
  
  
  § ACKNOWLEDGMENTS
   
  TY and JH were partially supported by U.S. National Science Foundation grant CMMI-2001262.  MD thanks the Velux Foundations for support under the Villum Experiment program (Project No. 00023059).
  We thank J. Gan for help with experiments, and M. Adda-Bedia and T. Jules for extensive discussions.  We also thank F. Lechenault, T. A. Witten, and J. Zhong for helpful discussions. 
  
  
  
  
  
  
  
  § NUMERICAL IMPLEMENTATION OF THE DEVELOPABLE MODEL
  
  
  The inextensible strip model is formulated as a two-point boundary value problem (BVP), which can be parametrically studied using the continuation package AUTO 07P <cit.>.
  This requires normalizing the length of the integral interval to unity.  This additional step is not explicitly shown in the following discussion, but can be achieved by replacing the s-derivatives (primes) with derivatives with respect to a normalized parameter s/s_max (if s ∈ [0,s_max]), thereby multiplying the right hand sides of all equations by s_max.  We retain the disk radius R in the descriptions below, although in our calculations we set it equal to unity for simplicity.
  
  To obtain inverted states, we begin with an annular sector of a flat disk subtending an angle ρ < π and bend this into one half of a conical frustum, the initial solution for numerical continuation.  The angle ρ is arbitrary and simply serves to create a nonflat conical starting point. Then the crease angle corresponding to a perfectly stiff crease is introduced by rotation of the ends of the sector, and ρ is increased to π. The crease stiffness is relaxed by replacing the boundary constraints on the crease angle with conditions on the moment.  Finally, the hole dimensions are adjusted.
  
  We employ Euler angles (ψ, θ, ϕ) to relate the director frame to a fixed Cartesian frame <cit.>,
  
  
      [ -N;  T;  B ]
      =
      [  cosϕ  sinϕ     0; -sinϕ  cosϕ     0;     0     0     1 ][  cosθ     0 -sinθ;     0     1     0;  sinθ     0  cosθ ][  cosψ  sinψ     0; -sinψ  cosψ     0;     0     0     1 ][ x; y; z ] .
  
  
  As κ_n=-N' ·T , κ_g=B' ·T , and τ_g=-B' ·N, we have
  
      κ_n=ϕ' +ψ' cosθ , 
      κ_g=  - θ' sinϕ + ψ' sinθcosϕ , 
      τ_g =θ' cosϕ +ψ' sinθsinϕ .
  
  Figure <ref> shows a Cartesian coordinate system (x,y,z) and a sequence of rotations applied to deform the annular sector into one conical half of a creased structure with axis z and bisected by the x-z plane.  The origin is at 12(r (0)+ r (ρ R) ).  
  
  
  
  For circular holes, the limit V of the generator coordinate v can be obtained explicitly in terms of the function η and the annular radii a and R <cit.>. 
  However, for elliptical holes, V has a complicated dependence on the backbone coordinate s and must be represented by an implicit function χ(V, s, η ;  a, b, R)=0.  Using η = β and λ=s/R, this is
  
      χ(s, η, V(s, η)) = V^2 +(R^2-2V R -b^2) sin^2 β
      
      	+ [(b/a)^2-1] [ Vcos(β + λ) + Rsinβsinλ]^2   ,
  
  which simplifies considerably for circular holes, for which a=b.
  
  
  
  
  
  
  We treat V as an independent variable and turn the algebraic constraint χ=0 into a differential equation.
  Denoting explicit partial derivatives by subscripts, we have V'= V_s + V_ηη'
  with V_s = - χ_s / χ_V and V_η = - χ_η / χ_V obtained through implicit differentiation <cit.>.
  We further differentiate the algebraic constitutive law in (<ref>) and combine with (<ref>) to obtain a first order ordinary differential equation (ODE) for κ_n and a second order ODE for η, and introduce another variable Ω   (=η') to convert the latter to two first order ODEs.
  The system is made autonomous by adding a trival differential equation s'=1.
  Constants “lost” through differentiation are added back using additional boundary conditions.
  All of this is combined with equations (<ref>-<ref>), (<ref>), and r'(s)=T to form the full system,
  
      F_1'-κ_n F_2 +κ_g F_3=0  , 
         F_2'+κ_n F_1 -κ_n η F_3=0  , 
         F_3'+κ_n η F_2 -κ_g F_1=0  ,  
         M_1'-κ_n M_2 +κ_g M_3=0  , 
         M_2'+κ_n M_1-κ_n η M_3 -F_3=0  , 
         M_3'+κ_n η M_2 -κ_g M_1+F_2=0  ,  
         η' = Ω , 
         (AE-C^2) Ω ' = (CB-AI) Ω+ AG -CJ  , 
         (AE-C^2) κ_n ' = (IC-BE) Ω + JE -GC  , 
         V'=-χ_s/χ_V-χ_η/χ_VΩ , 
         ψ' =(τ_g sinϕ + κ_g cosϕ)/ sinθ , 
         θ ' =τ_g cosϕ  - κ_g sinϕ , 
         ϕ ' =κ_n  - (κ_g cosϕ + τ_g sinϕ) / tanθ , 
         x'=-sinψcosϕ -cosψsinϕcosθ , 
         y'=cosψcosϕ -sinψsinϕcosθ , 
         z'=sinθsinϕ , 
         s'=1  ,
  
  in which 
  
      A   =Y_κ_n κ_n W  , 
      
       B   = Y_κ_n η W + Y_κ_n W_η + Y_κ_n W_V V_η , 
      
       C   = Y_κ_n η' W + Y_κ_n W_η' ,
      
      I   =Y_η' η W + Y_η' W_η +Y_η' W_V V_η
         +Y_η  W_η' +Y W_η' η +Y W_η' V V_η ,
      
      E   =Y_η' η' W + 2 Y_η' W_η' +W_η' η'Y  , 
      
       J   = η' M_1 -F_2 + κ_g (M_1 - η M_3) - Y_κ_n W_V V_s  ,
      
      G   =Y_η W + Y W_η +Y W_V V_η
         - κ_n M_1  - Y_η' W_V V_s - Y W_η' V V_s ,
  
  
  
  
  
  
  Figure <ref> shows the sequence of continuation steps we use to obtain inverted solutions.  
  A flat annular sector subtending an angle ρ < π with a small circular hole a ≪ 1 can be bent into the conical frustum of Figure <ref>b.  This starting solution for continuation is
  
      F_1=0 , F_2=0 , F_3=0 , 
         M_1=0 , M_2=lna/R , M_3=- √(π^2/ρ  ^2 - 1) lna/R , 
         κ_n=1/R√(π^2/ρ  ^2 - 1) , η=0 , η '=0 , 
         V= R -a   ,
         ψ= -π/2 + π/ρ Rs  , θ= -sin^-1ρ/π ,ϕ=0 , 
         x=ρ R/πsin (π/ρ Rs ) , y=-ρ R/πcos (π/ρ Rs ) ,z=0   ,
  
  with s ∈[0   , ρ R ].
  The ends of the sector are rotated about themselves to introduce the rest crease angle γ_0 (Figure <ref>c).  The boundary conditions are then 
  
      F_z (0)=0 , F_y (0)=0 , 
         M_x (0)=0 ,  
         κ_n (0) (1+η ^2 (0))^2/η' (0) + κ_g (1+η^2 (0)) W(η' (0) , η (0)) 
            - η(0) M_1 (0) - M_3 (0)=0 , 
         η(0)=0 , η(ρ R)=0 , 
         χ(V(0),0,η(0))=0  ,
         ψ(0)=-π/2 , ϕ (0)=- (π/2 - γ_0/2) ,
         ψ(ρ R)=π/2 , ϕ (ρ R)= (π/2  - γ_0/2) , 
         x(0) =0  ,z(0)=0,x(ρ R) =0 , z(ρ R)=0 , 
         y(0)+y(ρ R)=0 , 
         s(0)=0 ,
  
  where symmetry dictates the conditions on F_y(0)=F(0) ·y, F_z(0)=F(0) ·z, and M_x(0)=M(0) ·x.
  At this point the boundary conditions constrain the crease angle, so we have a perfectly stiff crease.
  In the next step, we increase ρ to π by a simple rescaling of s/s_max.
  Then the actual crease stiffness is introduced by replacing the two boundary conditions for ϕ above with moment conditions
  
      M_3 (0) =K_c R/D (1-b/R) sin[ π - γ_0 + 2 ϕ (0) ]   , 
         M_3 (π R) =K_c R/D (1-b/R) sin[ π - γ_0 - 2 ϕ (π R) ]   ,
  
  and decreasing K_c R/D from a large value to the real crease stiffness. 
  Finally, the hole dimensions a and b are continued to their real values.
  
  The geometry of the strip is reconstructed as 
  
  
      X(s,v)   =r(s)+v[B(s)+η(s)T(s)]  ,
         =(x-v [η (sinψcosϕ +cosψsinϕcosθ)-sinθcosψ]) x
          +(y+v [η (cosψcosϕ -sinψsinϕcosθ)+sinθsinψ ]) y
          +(z+v [ηsinθsinϕ +cosθ ]) z ,
  
  
  with v ∈ [0,V]. 
  
  The edge of regression is 
  
      c(s)=r(s)+ sinβ/β'-κ_gB(s)+η(s)T(s)/|B(s)+η(s)T(s)|
      
      = r(s)- B(s)+η(s)T(s)/η'+κ_g(1+η^2) .
  
  Its first derivative is c'(s)=ηκ_g^2  (1+η^2) + η”+3 κ_g ηη' [η'+κ_g (1+η^2)]^2 (B +ηT).
  At “cylindrical” points η'=-κ_g (1+η^2), the edge of regression goes off to infinity, and the mean curvature is constant along the local generator. 
  At “conical” points η”=-3 κ_g ηη' - ηκ_g ^2 (1+η^2), the edge of regression has a cusp.
  The generators can be mapped onto the flat annular sector with tangent t(s) and binormal b(s) using
  
      X(s,v)   =r(s)+v[b(s)+η(s)t(s)]  , 
          =(Rsinλ - v sinλ  + v ηcosλ )x
         + (- R cosλ + v ηsinλ + v cosλ) y ,
  
  and the edge of regression onto it using
  
      c(s)   =r(s)- b(s)+η(s)t(s)/η'+κ_g(1+η^2) , 
         = ( Rsinλ +sinλ  - ηcosλ/η'+κ_g(1+η^2))x
          -  (Rcosλ + ηsinλcosλ/η'+κ_g(1+η^2)) y .
  
  
  
  
  
  
  
  
  
  § FINITE ELEMENT SIMULATIONS
  
  
  Simulations were performed in the commercial finite element (FE) software COMSOL Multiphysics 5.4.  We used quadratic shell elements with a linear elastic Hookean material and geometrically nonlinear kinematic relations, and searched for solutions with the default stationary solver that implements the nonlinear Newton method. Mesh refinement studies were undertaken to ensure convergence of the results.  Symmetries of the disk were exploited so that only one quarter of the domain required simulation.
  Aspects of the simulations are illustrated schematically in Figure <ref>.
  
  
  
  The disk dimensions were set to R=75mm and t=127  μm, and the material properties used were a Young's modulus E=3.63 GPa and Poisson's ratio ν=0.4, consistent with reported material data and previous measurements on polyester sheets.
  
  The simulations consist of two steps.  First, a temperature-induced folding angle is set via a coupled phase-field-like model, in which the effect of a through-thickness temperature gradient is confined to within the crease region. Then, with the two outer crease ends vertically constrained, a downward force is applied at both inner crease ends to deform the structure from the folded configuration to the inverted configuration.
  
  The temperature-induced creasing follows the method introduced in <cit.>. 
  An effective thermal expansion coefficient ρ_T takes the value 1 K^-1 in a strip of width β t and 0 elsewhere.  In the present work, we use β=20, so that the width of the crease is comparable with that of the plastic region of real creases in mylar sheets <cit.>.
  The plate is subjected to a through-thickness temperature difference Δ T and responds with localized bending in the crease region.
  Empirically, we have found that for folded rectangular sheets, the rest crease angle γ_0 is given approximately by
  
      γ_0 = π - 1.4 ρ_T Δ T β  .
  
  By adjusting Δ T we are able to approximate a desired rest crease angle.  In the present work, while aiming for a γ_0 = 45^∘ we achieved  γ_0 = 48 ± 1^∘ with most of the deviation occurring near the edges of the crease. 
  
  During indentation, localized buckling can occur for some hole geometries.  To prevent this, foundation springs in the x direction were attached to the crease ends at the beginning of indentation, whose stiffnesses decrease linearly with the vertical indentation depth so that their effects vanish before snapping occurs.
  
  We estimate an effective crease stiffness for comparison with the developable model by ignoring boundary and ridge effects.  Let the moment per unit length of crease be
  
      M = ∫^t/2_t/2 z σ_xx dz  ,
  
  and the stress-strain relation be
  
      σ_xx = E/1-ν^2 (ϵ_xx-ϵ^0_xx)  ,
  
  where ϵ^0_xx is the “rest strain” due to thermal expansion.  The strains are related to the crease angle by
  
      ϵ_xx = π - γ/β t z  .
  
  Thus, the resultant moment on one half of the crease is
  
      M = - (R-b) D/β t (γ - γ_0)   .
  
  We can thus identify D/β t with K_c in the developable model, and R/β t with the dimensionless crease stiffness K_cR/D, where β = 20 in the present work.
  
  
  
  
  
  
  
  § DETAILS OF DEVELOPABLE SOLUTIONS
   
  
  Figures <ref> and <ref> provide details of the geometry of configurations from Figures <ref> and <ref>, respectively.
  
  
  
  
  
  
  
  
  
  
  
  unsrt
  
  10
  
  AmirbayatHearle86-1
  J. Amirbayat and J. W. S. Hearle.
  The complex buckling of flexible sheet materials— Part I.
    Theoretical approach.
  International Journal of Mechanical Sciences, 28(6):339–358,
    1986.
  
  AmirbayatHearle86-2
  J. Amirbayat and J. W. S. Hearle.
  The complex buckling of flexible sheet materials— Part II.
    Experimental study of three-fold buckling.
  International Journal of Mechanical Sciences, 28(6):359–370,
    1986.
  
  BenAmarPomeau97
  M. Ben Amar and Y. Pomeau.
  Crumpled paper.
  Proceedings of the Royal Society of London A, 453:729–755,
    1997.
  
  Chaieb98
  S. Chaïeb, F. Melo, and J.-C. Géminard.
  Experimental study of developable cones.
  Physical Review Letters, 80(11):2354–2357, 1998.
  
  CerdaMahadevan98
  E. Cerda and L. Mahadevan.
  Conical surfaces and crescent singularities in crumpled sheets.
  Physical Review Letters, 80(11):2358–2361, 1998.
  
  MoraBoudaoud02
  T. Mora and A. Boudaoud.
  Thin elastic plates: On the core of developable cones.
  Europhysics Letters, 59(1):41–47, 2002.
  
  lobkovsky1995scaling
  A. Lobkovsky, S. Gentges, H. Li, D. Morse, and T. A. Witten.
  Scaling properties of stretching ridges in a crumpled elastic sheet.
  Science, 270(5241):1482–1485, 1995.
  
  DiDonna02
  B. A. DiDonna.
  Scaling of the buckling transition of ridges in thin sheets.
  Physical Review E, 66:016601, 2002.
  
  LiangWitten05
  T. Liang and T. A. Witten.
  Crescent singularities in crumpled sheets.
  Physical Review E, 71:016612, 2005.
  
  FarmerCalladine05
  S. M. Farmer and C. R. Calladine.
  Geometry of “developable cones”.
  International Journal of Mechanical Sciences, 47:509–520,
    2005.
  
  Nasto13
  A. Nasto, A. Ajdari, A. Lazarus, A. Vaziri, and P. M. Reis.
  Localization of deformation in thin shells under indentation.
  Soft Matter, 9:6796–6803, 2013.
  
  ChopinKudrolli16
  J. Chopin and A. Kudrolli.
  Disclinations, e-cones, and their interactions in extensible sheets.
  Soft Matter, 12:4457–4462, 2016.
  
  Yang18
  Y. Yang, M. A. Dias, and D. P. Holmes.
  Multistable kirigami for tunable architected materials.
  Physical Review Materials, 2:110601(R), 2018.
  
  Moshe19
  M. Moshe, E. Esposito, S. Shankar, B. Bircan, I. Cohen, D. R. Nelson, and M. J.
    Bowick.
  Kirigami mechanics as stress relief by elastic charges.
  Physical Review Letters, 122:048001, 2019.
  
  Elder19
  T. Elder, D. Rozairo, and A. B. Croll.
  Origami inspired mechanics: Measuring modulus and force recovery with
    bent polymer films.
  Macromolecules, 52:690–699, 2019.
  
  Witten07
  T. A. Witten.
  Stress focusing in elastic sheets.
  Reviews of Modern Physics, 79:643–675, 2007.
  
  Das07
  M. Das, A. Vaziri, A. Kudrolli, and L. Mahadevan.
  Curvature condensation and bifurcation in an elastic shell.
  Physical Review Letters, 98:014301, 2007.
  
  Schroll11
  R. D. Schroll, E. Katifori, and B. Davidovitch.
  Elastic building blocks for confined sheets.
  Physical Review Letters, 106:074301, 2011.
  
  Mellado11
  P. Mellado, S. Cheng, and A. Concha.
  Mechanical response of a self-avoiding membrane: Fold collisions
    and the birth of conical singularities.
  Physical Review E, 83:036607, 2011.
  
  Balkcom09
  D. J. Balkcom, E. D. Demaine, M. L. Demaine, J. A. Ochsendorf, and Z. You.
  Folding paper shopping bags.
  In Origami4: Proceedings of the 4th International Meeting of
    Origami Science, Math, and Education (OSME 2006), pages 315–334, 2009.
  
  Tallinen10
  T. Tallinen, J. Ojajärvi, J. A. Åström, and J. Timonen.
  Scaling behavior in non-Hookean compression of thin-walled
    structures.
  Physical Review Letters, 105:066102, 2010.
  
  Guven13dipoles
  J. Guven, J. A. Hanna, O. Kahraman, and M. M. Müller.
  Dipoles in thin sheets.
  European Physical Journal E, 36:106, 2013.
  
  CambouMenon11
  A. D. Cambou and N. Menon.
  Three-dimensional structure of a sheet crumpled into a ball.
  Proceedings of the National Academy of Sciences,
    108(36):14741–14745, 2011.
  
  Thompson15
  J. M. T. Thompson.
  Advances in shell buckling: Theory and experiments.
  International Journal of Bifurcation and Chaos, 25(1):1530001,
    2015.
  
  Timounay20
  Y. Timounay, R. De, J. L. Stelzel, Z. S. Schrecengost, M. M. Ripp, and J. D.
    Paulsen.
  Crumples as a generic stress-focusing instability in confined sheets.
  Physical Review X, 10:021008, 2020.
  
  Hamm04
  E. Hamm, B. Roman, and F. Melo.
  Dynamics of developable cones under shear.
  Physical Review E, 70:026607, 2004.
  
  walsh2011weakening
  L. Walsh, R. Meza, and E. Hamm.
  Weakening of a thin shell structure by annihilating singularities.
  Journal of Physics D: Applied Physics, 44(23):232002, 2011.
  
  witten09spontaneous
  T. A. Witten.
  Spontaneous free-boundary structure in crumpled membranes.
  The Journal of Physical Chemistry B, 113(12):3738–3742, 2009.
  
  korte10triangular
  A. P. Korte, E. L. Starostin, and G. H. M. van der Heijden.
  Triangular buckling patterns of twisted inextensible strips.
  Proceedings of the Royal Society A, 467(2125):285–303, 2010.
  
  dias2012geometric
  M. A. Dias, L. H. Dudte, L. Mahadevan, and C. D. Santangelo.
  Geometric mechanics of curved crease origami.
  Physical Review Letters, 109(11):114301, 2012.
  
  Couturier13
  E. Couturier, J. Dumais, E. Cerda, and E. Katifori.
  Folding of an opened spherical shell.
  Soft Matter, 9:8359–8367, 2013.
  
  Seffen12
  K. A. Seffen.
  Compliant shell mechanisms.
  Philosophical Transactions of the Royal Society A,
    370:2010–2026, 2012.
  
  seffen18spherical
  K. A. Seffen.
  Spherical images and inextensible curved folding.
  Physical Review E, 97(2):023004, 2018.
  
  hanna2014waterbomb
  B. H. Hanna, J. M. Lund, R. J. Lang, S. P. Magleby, and L. L. Howell.
  Waterbomb base: a symmetric single-vertex bistable origami mechanism.
  Smart Materials and Structures, 23(9):094009, 2014.
  
  silverberg2014using
  J. L. Silverberg, A. A. Evans, L. McLeod, R. C. Hayward, T. Hull, C. D.
    Santangelo, and I. Cohen.
  Using origami design principles to fold reprogrammable mechanical
    metamaterials.
  Science, 345(6197):647–650, 2014.
  
  Waitukaitis15
  S. Waitukaitis, R. Menaut, B. G. Chen, and M. van Hecke.
  Origami multistability: From single vertices to metasheets.
  Physical Review Letters, 114:055503, 2015.
  
  lechenault15generic
  F. Lechenault and M. Adda-Bedia.
  Generic bistability in creased conical surfaces.
  Physical Review Letters, 115(23):235501, 2015.
  
  andrade2019foldable
  I. Andrade-Silva, M. Adda-Bedia, and M. A. Dias.
  Foldable cones as a framework for nonrigid origami.
  Physical Review E, 100(3):033003, 2019.
  
  gillman2018truss
  A. Gillman, K. Fuchi, and P. R. Buskohl.
  Truss-based nonlinear mechanical analysis for origami structures
    exhibiting bifurcation and limit point instabilities.
  International Journal of Solids and Structures, 147:80–93,
    2018.
  
  walker18shape
  M. G. Walker and K. A. Seffen.
  On the shape of bistable creased strips.
  Thin-Walled Structures, 124:538–545, 2018.
  
  Walker20
  M. G. Walker.
  Mechanics of generically creased disks.
  Physical Review E, 101:043001, 2020.
  
  lechenault2014mechanical
  F. Lechenault, B. Thiria, and M. Adda-Bedia.
  Mechanical response of a creased sheet.
  Physical Review Letters, 112(24):244301, 2014.
  
  francis2013origami
  K. C. Francis, J. E. Blanch, S. P. Magleby, and L. L. Howell.
  Origami-like creases in sheet materials for compliant mechanism
    design.
  Mechanical Sciences, 4(2):371–380, 2013.
  
  starostin2007shape
  E. L. Starostin and G. H. M. van der Heijden.
  The shape of a Möbius strip.
  Nature Materials, 6(8):563–567, 2007.
  
  starostin2015equilibrium
  E. L. Starostin and G. H. M. van der Heijden.
  Equilibrium shapes with stress localisation for inextensible elastic
    Möbius and other strips.
  Journal of Elasticity, 119(1-2):67–112, 2015.
  
  MooreHealey18
  A. Moore and T. Healey.
  Computation of elastic equilibria of complete Möbius bands
    and their stability.
  Mathematics and Mechanics of Solids, 24(4):939–967, 2018.
  
  yu2019bifurcations
  T. Yu and J. A. Hanna.
  Bifurcations of buckled, clamped anisotropic rods and thin bands
    under lateral end translations.
  Journal of the Mechanics and Physics of Solids, 122:657–685,
    2019.
  
  badger2019normalized
  J. C. Badger, T. G. Nelson, R. J. Lang, D. M. Halverson, and L. L. Howell.
  Normalized coordinate equations and an energy method for predicting
    natural curved-fold configurations.
  Journal of Applied Mechanics, 86(7):071006, 2019.
  
  todres2015translation
  R. E. Todres.
  Translation of W. Wunderlich’s “On a developable
    Möbius band”.
  Journal of Elasticity, 119(1-2):23–34, 2015.
  
  dias2015wunderlich
  M. A. Dias and B. Audoly.
  “Wunderlich, meet Kirchhoff”: A general and unified description
    of elastic ribbons and thin rods.
  Journal of Elasticity, 119(1-2):49–66, 2015.
  
  doedel2007auto
  E. J. Doedel, R. C. Paffenroth, A. R. Champneys, T. F. Fairgrieve, Y. A.
    Kuznetsov, B. E. Oldeman, B. Sandstede, and X. Wang.
  AUTO-07P: Continuation and bifurcation software for ordinary
    differential equations.
  
   , 2007.
  
  love1927treatise
  A. E. H. Love.
  A treatise on the mathematical theory of elasticity.
  Dover, New York, 1927.
  
  jules19local
  T. Jules, F. Lechenault, and M. Adda-Bedia.
  Local mechanical description of an elastic fold.
  Soft Matter, 15(7):1619–1626, 2019.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  yu2021bistability
  T. Yu.
  Bistability and equilibria of creased annular sheets and strips.
  arXiv:2104.09704, 2021.
  
  
  
  