## Abstract

Self-folding systems, which can transform autonomously from a flat sheet into a 3D machine, provide opportunities for rapidly fabricable robots that are deployable on-demand. Existing self-folding fabrication processes convert fold patterns into laminated structures that respond to external stimuli, most commonly heat. However, demonstrations of these approaches have been generally limited to simple fold patterns with little ambiguity in folding configuration, and the reliability of self-folding drops drastically with the fold pattern complexity. In this paper, we explore methods of biasing a symmetric fold pattern, the origami hyperbolic paraboloid (hypar), to fold into one of the two possible configurations. The biasing methods are simulated using a bar-and-hinge inspired self-folding model that defines a single fold as a bending beam and the hypar crease pattern as an elastic spring network. Simulation results are also verified on physical samples. Based on these results, three techniques to bias the hypar by manipulating the target fold angles are proposed and tested. The results show that biasing a self-folding pattern can increase folding accuracy from 50% (purely random) to 70% and provide insights for improving the reliability of future self-folding systems with complex fold patterns.

## 1 Introduction

Self-folding is a self-assembly process by which a structure can transform autonomously from a flat sheet into a 3D shape [1–4]. Inspired by origami principles, self-folding presents several benefits over traditional methods of rapid fabrication for 3D objects, including less material waste and improved cost-efficiency [4]. These structures are ideal for robotic systems since multiple materials, including actuators and sensors, can be integrated directly into the sheet [5], and what is essentially a flexible printed circuit board can then perform robotic tasks with little required wiring and manual assembly. This logistical advantage has increased interest in self-folding systems, as the devices can be transported in the flat state and deployed autonomously into complex geometries at millimeter to centimeter scales [6–10]. Recent studies have explored different actuation approaches to self-folding, including hydrogels [11], infrared light absorption [12], pH [13], magnetic fields [14,15], capillary forces [16], and shape-memory alloys [17]. In general, these approaches place active and passive materials at the folds and faces, respectively. When an external signal is applied, the active material actuates along the folds, applying rotational transformations to the faces and folding the structure into its 3D form.

A significant challenge in all of these approaches is consistency and reliability of folding. While most self-folding techniques have been demonstrated reliably on simple patterns with a small number of folds, the error rate grows drastically when the fold pattern increases in complexity. The main problem lies in the pattern’s flat, singular configuration. From this state, multiple folding pathways branch out in trajectories that the fold pattern can take into its various final folded states. Theoretical analyses [18,19] indicate that the number of branches (i.e., the number of different configurations that a pattern is able to fold into) increases exponentially with the number of folds, and even indicating the direction of folding through a mountain-valley assignment for every single fold in the pattern is not sufficient to differentiate a particular target folded state from the other possibilities [20]. Thus, for the purposes of remote deployment, the question of how to best program a fold pattern to self-fold into the target shape consistently remains a problem of great interest [21].

Existing approaches to fold pattern biasing fall into two main categories: biasing the folds and biasing the vertices. For biasing the folds, Tachi and Hull [20] first analyzed the configuration spaces of various fold patterns and observed that the forces applied at each of the folds could be manipulated to make the pattern prefer one branch over another. In particular, by choosing a configuration space force whose dot product with the tangent to the desired branch is positive and that is normal to all other branches, a single branch can be chosen. The force is then implemented by mapping the individual force values to the folds on the physical pattern. This approach was demonstrated successfully on a simple degree-6 vertex. Unfortunately, the authors also showed that there exist some patterns for which finding this ideal driving force is not possible. Furthermore, as the complexity of the pattern increases, mapping the entire configuration space of the pattern quickly becomes intractable. Thus, Stern et al. [19,22] extended this idea by proposing a heuristic whereby a force tangent to the desired branch is used. This approach significantly reduces the complexity of computing the driving force and further has the advantage that only a subset of the folds may need to be actuated at all [22]. The approach was demonstrated successfully in simulation, although the authors did note that the success rate decreased for fast systems with inertia. To improve self-folding consistency, biasing the vertices can be used. In particular, Chen and Santangelo [18] conducted numerical analyses of single-vertex patterns and discovered that in the large majority of cases, the mountain-valley assignment together with the orientation of the center vertex (popped up or down) was sufficient to uniquely identify the target folded state. Kang et al. [23] verified this strategy experimentally on a degree-6 vertex and a bird pattern.

These approaches are designed around rigid origami. In particular, the configuration space analysis and resulting design approaches rely on an assumption that kinematic constraints are maintained and that the only potential energy in the system comes from the fold bending deformations. Practically, because self-folding structures are fabricated from thin films, faces can deform significantly away from the rigid assumption, and there exist, indeed, many examples of non-rigid origami patterns that may be useful in deployable systems.

To address this gap in understanding, we study the self-foldability of a non-rigid origami pattern, the origami hyperbolic paraboloid (hypar) [24,25]. The hypar is especially of interest because of its symmetry and bistability. The hypar is folded from concentric squares and their diagonals, as shown in Fig. 1(a). When folded, the hypar pattern’s discrete pleats converge to a smooth global saddle shape [26,27], modeling a negative Gaussian curvature surface with hyperbolic and parabolic cross sections [28] (Fig. 1(b)). Because it is not rigidly foldable, the hypar achieves its three-dimensional shape through elastic bending deformations not explicit in the fold pattern; these hidden bending degrees-of-freedom also give rise to the hypar’s symmetrically bistable behavior [29]. Interestingly, the hypar’s final configuration is not unique defined by its mountain-valley assignment. The hypar has two symmetric stable folded states: in one case, one pair of diagonals is pointed up while the other is pointed down; in the other case, the opposite is true. This symmetrically bistable characteristic makes programing the hypar to self-fold into a particular configuration challenging.

Our work explores self-folding of the origami hypar in both simulation and physical experiments. The main contributions of this work include:

An extension of existing bar-and-hinge models to self-folding structures.

Simulation and analysis of simple patterns (including a single-fold pattern and a degree-4 vertex) and the hypar.

Simulated and experimental exploration of fabrication strategies for biasing the hypar to fold into one stable state or another.

## 2 Self-Folding Fabrication Overview

Most self-folding fabrication processes take a multi-material approach, where the folds of the fold pattern are fabricated from some active material, and the faces are fabricated from a semi-rigid, passive material. In this work, we model the laminate-based approach used in Refs. [7,31–33], wherein three layers of plastic film are adhered together to form the folds and faces, as outlined in Fig. 1, although the model is extensible to other self-folding methods. In particular, a thermoresponsive polyvinyl chloride (PVC) film is sandwiched between two structural layers made of polyethylene terephthalate (PET). Gaps cut into the structural layer expose the PVC such that the entire structure folds along these gaps when the structure is heated. This deformation is preserved as the laminate cools, resulting in a permanently folded three-dimensional structure.

The self-folding structure is designed and fabricated using a multi-step process. The original fold pattern (Fig. 1(a)) is first converted into two mirror-image cut patterns forming the structural layers (Fig. 1(c)). Gaps corresponding to valley folds are placed on the left image, and gaps corresponding to mountain folds are placed on the right image. Because the size of the gap dictates the final fold angle (refer Sec. 3.1), precise alignment of all the faces must be maintained during fabrication. To do this, we first adhere PET film (Grafix Dura-Lar.002-Ultra-Clear Film) onto a sacrificial layer of painter’s tape (Fig. 1(d)). The PET film is also covered with silicone glue (McMaster part 7615A41, 0.051 mm thick). We then cut the gap pattern onto the layers using a laser cutter (Universal Laser Systems PLS 4.75), cutting through the PET film but not through the painter’s tape. The gaps are removed from the painter’s tape, leaving only the faces behind (Fig. 1(e)). Next, 0.025 mm thick PVC film is laid on the left side of the gap pattern and then entire structure is folded across the centerline, sandwiching the PVC between the two structural layers and aligning all of the faces (Fig. 1(g)). Finally, the painter’s tape is removed, leaving behind the self-folding structure. The resulting fold pattern is 0.229 mm thick. For self-folding, the structure is then submerged into near-boiling water (90 °C), where it folds into shape within 500 ms (Fig. 1(h)). Figure 1(f) shows a cross-sectional view of the laminate structure before and after folding.

## 3 Bar-and-Hinge Model for Self-Folding Structures

To model this self-folded structure, we adapt the nonlinear elastic bar-and-hinge method [34–36]. In this model, a crease pattern is modeled as an elastic spring network composed of bar elements, folding hinges, and bending hinges. Each face consists of bar elements along each of its edges. These elements are nonlinear springs, which respond to the stretching of the face. Additional bars across the diagonals of the face track the shear deformation. The placement and connectivity of diagonal bars is the subject of multiple investigations. We use the N5B8 model from Ref. [35]. Bending of the structure is tracked as torsional springs positioned at the folds and face diagonals (called “bends”). Note that the stiffness of the folds may differ from the face diagonals to account for creasing, material fatigue, or, in the case of self-folded structures, gaps cut into the structural layers.

**u**be a vector of nodal displacements for each of the nodes in the spring representation of the pattern. Then,

*U*

_{s}(

**u**) is the strain potential energy stored in the bar elements,

*U*

_{f}(

**u**) is the potential energy stored in the folds, and

*U*

_{b}(

**u**) is the potential energy stored in the bends. Then the total potential energy stored in the pattern is

### 3.1 Fold Model for Self-Folding Structures.

We extend upon the standard bar-and-hinge model by incorporating a new model at the folds that is based on the self-folding fabrication process. In particular, when the self-folding structure is exposed to heat, the PVC in the middle layer shrinks, while the PET at the outer layers does not. This has the net effect of changing the equilibrium angle (and thus the potential energy profile) of the folds. If the new potential energy profile is known, then we can simulate the self-folded shape by iteratively computing the equilibrium configuration of the fold pattern as the PVC shrinks. These equilibria configurations are solved for using a damped Newton–Raphson algorithm. Since we are particularly interested in challenging (i.e., symmetric) fold patterns, and these patterns have a bifurcation at when the pattern is flat, the simulation may be unstable when the fold angles are close to 0. An adaptive gap width increment, whereby gap width is increased by a smaller amount when the Newton–Raphson solver fails to converge and by a larger amount when it succeeds, is used to improve stability of the solver near this bifurcation.

Then it simply remains to model a single fold. We model a single fold as a bending beam constrained by the thickness of the rigid material. In particular, consider a laminate structure with structural PET layers of thickness *t*_{s} and a shrinking PVC layer of thickness *t*_{g}. Gaps are cut into one layer of PET so that the PVC is exposed by a gap width of *g*. Now, when the structure folds, the PVC shrinks by a factor *λ*, resulting in the sheet of material *λg* wide that is supposed to cross a *g*-width gap. Since the uncut PET layer remains, the structure bends.

To solve for the resulting bend angle and stiffness of the fold, we use the following assumptions:

The folds are independent: bending of one fold does not affect bending of another. This approximation is reasonable for small gap widths. For large gap widths, however, a fold may see some small influence from neighboring folds since material is shared at the vertex.

The fold is 1D: bending of a fold is constant over the entire fold. This is true for rigidly foldable models but is not true for models like the hypar. To address this issue, we can segment a pattern into shorter folds such that each individual small fold can be approximated as having a constant fold angle.

The gap is symmetric: thus, it is only necessary to model half of the fold. Furthermore, the fold can be modeled as a cantilever beam fixed at one end.

The structural layers are rigid, and there is no delamination at the faces.

We consider a bending beam of length *λg*/2 that is fixed at one end at an angle 0 and at the other end to a slanted, rigid plate. There are two cases to consider (Ref. Fig. 2), as determined by whether or not the rigid layers on the interior of the fold collide.

*Small gap sizes*(Fig. 2(a)): when the gap size is small, the final fold angle is dominated by the collision of the PET faces on the interior of the fold. In this case, for a fold with gap width*g*and material thickness*t*_{s}and*t*_{g}, the resulting fold angle is simply(2)$\theta =2tan\u22121(g/2ts+tg)$*Large gap sizes*(Fig. 2(b)): when the gap size is large, the fold angle is controlled by the curvature of the PVC layer, and the PET layers on the interior of the fold do not come into contact. In fact, there is an upper bound on the maximum fold angle for the laminate structure. Since fold angles can be quite high, we use a nonlinear beam bending model to analyze this case. Let*x*(*s*) and*y*(*s*) be the (*x*,*y*) coordinates of a point*s*distance along the length of the gap, where*s*∈ [0,*g*/2]. The curvature of the beam is governed by the moment equationWe assume that the fold is fixed such that(3)$d2ydx2=M(y)EI[1+(dydx)2]3/2$*dy*/*dx*= 0 at the center of the fold, and that the only attachment between the PVC and the PET faces is at the edge of the gap. That is, the only force being applied to the PVC fold is a point force at the tip. Assuming that this force is applied normally to the PET layer (i.e., that there is no shear between layers), the moment isIn addition, the system must satisfy the constraint that the uncut PET layer remains straight and of length(4)$M(y)=Py((y(L)\u2212y)+dydx(x(L)\u2212x))$*g*/2, that is,And, finally, the total length of the beam must be(5)$dydx|s=L=(g/2)2\u2212x2x$*λg*/2.The system can be solved using a numerical solver (in our case, matlab’s fsolve). Note that because the system is fully constrained, the actual values of

*P*,*E*, and*I*do not affect the final fold angle, although they do affect the stiffness of the fold. The final fold angle (which is simply twice the angle of the tip) is constant regardless of gap width. We will call this value*θ*_{max}.*Determining the case*: which case a fold falls into is then determined by taking the minimum of the computed fold angles. Thus the transition between the small gap case and the large gap case occurs at(6)$\theta =2tan\u22121(g/2ts+tg)=\theta max$*Fold stiffnesses*. The fold stiffnesses follow naturally from this model as the derivative of the force*P*applied at the edge of the gap with respect to the fold angle*θ*and can also be estimated numerically. Note that the stiffness is derived directly from the bending beam model and thus is the same for both the small and large gaps. Theoretically, forces should increase towards infinity in the direction of decreasing*θ*for small gaps since this configuration corresponds to self-collision. We ignore this issue in our simulation, although self-collision forces could be incorporated using barrier functions in a manner similar to Ref. [39].

### 3.2 Experimental Verification.

To check the validity of our assumptions and of the model, we conducted physical tests following the procedure in Sec. 2 and measured the resulting fold angles over a range of fold patterns and gap widths.

#### 3.2.1 Material Characterization.

The resulting material properties used are shown in Table 1. Using the individually measured Young’s modulus values, we computed the net Young’s modulus of the laminate as the weighted average of the individual layers, scaled by thickness. The total thickness of the laminate was computed as the sum of thicknesses of all the material layers (including the adhesive).

Sample | Thickness | Young’s modulus |
---|---|---|

PET | 0.051 mm | 2.9 GPa |

Unshrunk PVC | 0.025 mm | 3.8 GPa |

Shrunk PVC | 0.060 mm | 1.5 GPa |

Unshrunk laminate | 0.229 mm | 3.0 GPa |

Shrunk laminate | 0.264 mm | 2.4 GPa |

Sample | Thickness | Young’s modulus |
---|---|---|

PET | 0.051 mm | 2.9 GPa |

Unshrunk PVC | 0.025 mm | 3.8 GPa |

Shrunk PVC | 0.060 mm | 1.5 GPa |

Unshrunk laminate | 0.229 mm | 3.0 GPa |

Shrunk laminate | 0.264 mm | 2.4 GPa |

*t*

_{s}= 0.11 mm (i.e., the combined thickness of the PET and adhesive layers) and

*t*

_{g}= 0.06 mm. Thus, the theoretical transition between the small gap and large gap cases is

*g*= 0.81 mm, with a corresponding

*θ*

_{max}= 2.36 rad.

#### 3.2.2 Single-Fold Model.

To verify the basic fold angle model, we fabricated and measured physical samples according to the fabrication process outlined in Sec. 2. Two sets of samples were tested. In one set, we tested a simple bookfold pattern, which consists of two rectangular faces separated by a single fold. In the second set, we tested an accordion fold pattern, where four rectangular faces are separated by three folds of alternating mountain-valley assignment, in order to determine whether folds influence each other. For each set, we fabricated samples of varying face size and gap width. In particular, we tested rectangles 12.5 mm in length and 3.125 mm, 6.25 mm, or 12.5 mm in width. In this case, we use *length* to denote the dimension parallel to the fold, so all tested folds were also 12.5 mm long. Tested gap widths varied between 0 mm and 3 mm in increments of 0.25 mm. Figure 4 shows some of the self-folding accordion samples used.

The resulting fold angles were measured by imaging the self-folded patterns on their side and locating the vertices. Since many of the fold patterns had some twist (i.e., fold angles were different on side of a fold versus the other), this process was repeated for both sides of the fold pattern and averaged into a single-fold angle value. Figure 5 shows the results. The individual error bars show the mean, maximum, and minimum angle measured for each pattern. Since there was no observable difference between the fold angles on the accordion folds pattern and fold angles on the bookfold pattern, the two data sets were combined for analysis. The black line shows the results of our model. The portion on the left is the small angle model, where the interior layers of PET collide, while the portion on the right is the constant maximum fold angle of 2.36 rad. As predicted by our model, the physical samples clearly show a trend of increasing fold angle with increasing gap width and level out with a maximum angle of about 2.36 rad. Variance in the fold angles of the physical samples is likely due to fabrication error, in particular, slight misalignment between layers.

We can also observe horizontal scaling in the physical data compared to the model. In particular, physical samples tend to achieve a lower fold angle than predicted, and this discrepancy is larger for the wider faces than for the thinner faces. We expect that this is a result of delamination in the layers. Since the PET layers are attached to the PVC with a thin layer of glue, it is possible that the layers shift in some places, with the net effect of increasing the amount of PVC material that participates in the fold, effectively increasing the shrink ratio *λ* and reducing the fold angle. Using a least-squares fit, we computed that the shrink ratio is approximately 1.31*λ* for the 3.125 mm face width, 2.20*λ* for the 6.25 mm face width, and 2.32*λ* for the 12.5 mm face width. The different lines in Fig. 5 show the results with this scaling. Practically, this value must be measured experimentally, although these results do seem to indicate that the effect of face size is reduced when the face is large.

#### 3.2.3 Degree-4 Vertex.

As a simple extension of the basic self-folding model, we compare our model and simulation to physical prototypes of a degree-4 vertex fold (Fig. 6). This pattern is rigidly foldable, and the theoretical fold angles are well known. We fabricated degree-4 fold patterns from 25 mm × 25 mm squares with gap widths of both 0.5 mm and 1.25 mm. The gap widths were converted directly into fold angles using the scaled model in Sec. 3.1, yielding target fold angles of ±0.51 rad and ±1.15 rad, respectively, on all of the folds. Note that because of its kinematic constraints, a degree-4 vertex does not fold with all of the fold angles achieving the same value. We can therefore use this test as a check to ensure that relative fold stiffnesses are accurate and that the model accurately captures the relationships between folds.

We measured the resulting fold angles on each of the four folds of this pattern by tracking the edges of the fold pattern, similarly to Sec. 3.2.2. The edges used to measure angles are indicated in red lines shown in Fig. 6. Note that these measured angles are not exactly the same as the fold angles since they are not taken along a line perpendicular to the fold, but the measurements can still be used to verify the accuracy of the simulation since the angles were measured along the same edges. Three sample self-folding patterns were fabricated for each gap width, and the measured angles averaged across the three samples. The results are shown in Table 2. Both the physical and simulated patterns achieve the same mountain-valley assignment for each fold, and the fold angles of the physical samples match the simulated patterns to within 13%, with larger errors occurring for the 5 mm gap than the 1.25 mm gap. This result reflects the results of Fig. 5, which also showed larger variance in fold angle for the 0.5 mm gap than the 1.25 mm gap.

Gap (mm) | Measured angle (rad) | ||||
---|---|---|---|---|---|

θ_{1} | θ_{2} | θ_{3} | θ_{4} | ||

Sim | 0.5 | 2.56 | 1.53 | 2.56 | −1.52 |

1.25 | 1.83 | 1.35 | 1.84 | −1.34 | |

Exp | 0.5 | 2.23 | 1.46 | 2.19 | −1.45 |

1.25 | 1.85 | 1.39 | 1.98 | −1.42 |

Gap (mm) | Measured angle (rad) | ||||
---|---|---|---|---|---|

θ_{1} | θ_{2} | θ_{3} | θ_{4} | ||

Sim | 0.5 | 2.56 | 1.53 | 2.56 | −1.52 |

1.25 | 1.83 | 1.35 | 1.84 | −1.34 | |

Exp | 0.5 | 2.23 | 1.46 | 2.19 | −1.45 |

1.25 | 1.85 | 1.39 | 1.98 | −1.42 |

Notes: (Rows 1–2) Simulation results. (Rows 3–4) Experimental results.

#### 3.2.4 Hyperbolid Paraboloid (Hypar).

Finally, we tested a three-layer hyperbolic paraboloid (hypar) pattern. The same material parameters were used as in Sec. 3.2.3. Since the hypar is not rigidly foldable, the pattern was discretized along the pleat folds, with each of the outermost pleats having 30 segments. The pattern is shown in Fig. 7.

For ease of discussion, we label the vertices on the hypar crease pattern as follows: the corners of the hypar pattern are labeled 1–4 and written in subscripts. The pleats are numbered 0–3 from outside to inside and written in superscripts. Thus, the vertex $v10$ is the vertex in on the boundary in the top left corner, as shown in Fig. 7. The vertex $v43$ would be the vertex on bottom left corner of the innermost square. Then folds are denoted by hyphens. For example, $v10\u2212v11$ is the diagonal fold in the top left corner of the pattern.

The hypar pattern is symmetric and can fold into two saddle configurations. In each configuration, one set of diagonals moves up relative to the pattern and the other set of diagonals moves down. We call the configuration where the 1 and 3 corners are up the 1–3 configuration, and other the 2–4 configuration.

The results of simulating this pattern are shown in Fig. 7(b). Qualitatively, the simulation is able to predict the shape of the folded hypar. Furthermore, the simulation was equally able to fold into the 1–3 configuration and the 2–4 configuration depending on how the numerical solver within the simulation was initialized. In addition, the figures show the deviation of each fold from the programed fold angle. The majority of fold edges have deviations close to 0 rad, indicating that their actual fold angles are very close to the prescribed fold angles; quantitatively analyzing the final fold angles reveals that most errors lie between ±0.2 rad, with a 9.59% average angle difference.

*Parameter Randomization*. Our main goal in this work is to understand the reliability of self-folding and how to bias the hypar to fold into one configuration or another. We therefore performed an additional test where we introduced randomness into the pattern to simulate the effects of fabrication error. We are particularly interested in two sources of noise:

Vertex offset: During the fabrication of the self-folding laminate, misalignment of the two structural PET layers frequently occurs when the PVC film is inserted, resulting in unintentionally displaced nodes. To account for this misalignment error, we add small random error displacements to each node following a normal distribution $N(0,g/2)$ mm, where

*g*is the gap width. A check ensures that the absolute error does not exceed half the gap width.Fold angle variation: Misalignment errors not only affect node coordinates but also the equilibrium fold angle. Misalignment changes the area of the exposed PVC film, causing the fold angle to change as shown in Fig. 5. Thus, we adjust each fold angle target in the fold pattern by adding random errors following a normal distribution $N(0,\sigma a)$ rad, where

*σ*_{a}is the standard deviation of fold angles computed from our single-fold experiments. A check ensures that generated errors are positive and the target fold angle does not exceed*π*.

We simulated hypar patterns with a gap width of 0.5 mm. The pattern was simulated over 50 trials with the noise parameters shown in Table 3. The notation used is that *N*(*μ*_{n}, *σ*_{n}) indicates vertex offset where each vertex was displaced according to a normal distribution with mean *μ*_{n} and standard deviation *σ*_{n} in the *x* and *y* directions. A randomization setting of A(*μ*_{a}, *σ*_{a}) indicates that fold angle noise was introduced, where every angle was offset according to a normal distribution with mean *μ*_{a} and standard deviation *σ*_{a}. For a 0.5 mm gap width, we used *N*(0, 0.25) mm and *A*(0, 0.09) rad, which were found experimentally during our single-fold tests.

Biasing | Randomization parameters | 1–3 | 2–4 | No distinct configuration | p-Value |
---|---|---|---|---|---|

Unbiased | N(0, 0.25), A(0, 0.09) | 23 | 21 | 6 | 0.44 |

Center (1–3) | N(0, 0.25), A(0, 0.09) | 42 | 0 | 8 | <0.0001 |

Center (2–4) | N(0, 0.25), A(0, 0.09) | 0 | 45 | 5 | <0.0001 |

Pleat (1–3) | N(0, 0.25), A(0, 0.09) | 41 | 1 | 8 | <0.0001 |

Pleat (2–4) | N(0, 0.25), A(0, 0.09) | 0 | 46 | 4 | <0.0001 |

Diagonals (1–3) | N(0, 0.25), A(0, 0.09) | 50 | 0 | 0 | <0.0001 |

Diagonals (2–4) | N(0, 0.25), A(0, 0.09) | 0 | 50 | 0 | <0.0001 |

Biasing | Randomization parameters | 1–3 | 2–4 | No distinct configuration | p-Value |
---|---|---|---|---|---|

Unbiased | N(0, 0.25), A(0, 0.09) | 23 | 21 | 6 | 0.44 |

Center (1–3) | N(0, 0.25), A(0, 0.09) | 42 | 0 | 8 | <0.0001 |

Center (2–4) | N(0, 0.25), A(0, 0.09) | 0 | 45 | 5 | <0.0001 |

Pleat (1–3) | N(0, 0.25), A(0, 0.09) | 41 | 1 | 8 | <0.0001 |

Pleat (2–4) | N(0, 0.25), A(0, 0.09) | 0 | 46 | 4 | <0.0001 |

Diagonals (1–3) | N(0, 0.25), A(0, 0.09) | 50 | 0 | 0 | <0.0001 |

Diagonals (2–4) | N(0, 0.25), A(0, 0.09) | 0 | 50 | 0 | <0.0001 |

Notes: All biasing methods succeeded in changing the distribution of folded hypar configurations in simulation, with *p*-values less than 0.0001. Comparing the raw results, the diagonals biasing method produced the most consistent results with the fewest number of non-distinct foldings.

The effect of random noise on the hypar simulation is shown in Table 3. Both configurations were ultimately achieved, with no significant bias towards one configuration or another (*p*-value 0.44 as compared to a distribution with 50% chance of either configuration). This lack of bias confirms the symmetry of the pattern, demonstrating that the hypar is not intrinsically prone towards either configuration. As a comparison, 10 physical samples were also constructed and self-folded (Table 4). The result was 5 of each configuration, again confirming the symmetry of the pattern.

Biasing | Correct | Incorrect | p-Value |
---|---|---|---|

Unbiased | 5 | 5 | 0.5 |

Center | 11 | 10 | 0.5 |

Pleat | 14 | 6 | 0.059 |

Diagonals | 14 | 6 | 0.059 |

Biasing | Correct | Incorrect | p-Value |
---|---|---|---|

Unbiased | 5 | 5 | 0.5 |

Center | 11 | 10 | 0.5 |

Pleat | 14 | 6 | 0.059 |

Diagonals | 14 | 6 | 0.059 |

Note: While none of the results were statistically significantly different from the unbiased hypar, the pleat and diagonals biasing methods do show a greater number of correct configurations than the center and unbiased methods.

In up to 10% of the simulation trials, non-hypar final configurations were observed; the geometry of these incongruous trials were typically a flat, unfolded sheet or a crumpled form. These incongruities were consistent with empirical observations, in which the hypars did not self-fold perfectly in every single trial, especially those with severe misalignment errors. Because the hypar was able to achieve both configurations as well as some non-hypar configurations, the results corroborate the significance of random noise in the self-folding of the hypar: the final configuration is extremely sensitive to small nodal, angle, and force displacements, and some biasing mechanism is needed to produce reliable self-folding.

## 4 Biasing Mechanisms

In order to produce more reliable self-folding and deployment of origami structures, we investigate several biasing methods for the origami hypar, two that we designed specifically for the hypar and one inspired from existing theory [19]. Figure 8 shows a summary of the approaches, which are explained in the following subsections. Each row corresponds to a biasing method. For each method, the fold pattern design with target fold angles is shown in column 1. In column 2, the fabrication plan with the fold angles mapped to gap widths is shown. Columns 3 and 4 show the fabricated flat sheet and folded structure, respectively.

The simulation results for all biasing methods are summarized in Table 3, and the results of physical tests are in Table 4. For each randomization parameter and biasing method, 50 simulation trials and 10 physical tests of each configuration were run. For each result, we performed a one-tailed *z*-test to determine statistical significance, comparing against the hypothesized result of 50% for each configuration.

### 4.1 Adding and Removing Folds (Center Biasing).

In the first biasing method, we observed that the unbiased hypar pattern has two creases ($v13\u2212v33$ and $v23\u2212v43$) meeting at 90 deg in a cross pattern at the center. For this arrangement of creases, only one crease can fold at a time. Each of the hypar’s two configurations has a unique center diagonal fold assignment. In the 1–3 configuration, the $v23\u2212v43$ fold is a mountain fold, while the $v13\u2212v33$ fold remains flat; in the 2–4 configuration, the opposite pattern is observed, where the $v13\u2212v33$ fold is a mountain fold. Building upon these observations, we hypothesized that removing the flat folds and leaving only one fold in the center of the pattern could bias the final configuration of the hypar into whichever configuration used the fold that remained. The original unbiased crease pattern was thus modified to leave only one center diagonal, resulting in the crease pattern shown in Fig. 8 (row 2).

The results for center biasing indicate that while this strategy is highly successful in simulation, it fails to work in practice. In simulation, the hypars folded correctly for 84% (1–3) and 90% (2–4) of the trials; with a *p*-value of <0.0001 for both configurations, the results demonstrate that pleat biasing significantly biases the hypar. However, physical samples showed that 11/21 samples folded into the correct configuration, while 10/21 samples folded into the incorrect configuration. This distribution is essentially 50–50, showing that the center biasing had very little effect at all. This result is not unexpected. Only the center of the pattern was modified, but the outer edges of the pattern experience the most significant deformation during folding and the greatest resistance forces from the surrounding water bath. It thus makes sense that this biasing method should have minimal effect. This result indicates that additional considerations, for example, external forces, should be added to the simulation in future work.

### 4.2 Changing Fold Angle Along a Fold (Pleat Biasing).

In the second biasing method, we observed instead the changes in fold angle that occur along a fold of the hypar. Since the hypar is not a rigidly foldable pattern, the panels bend and twist, and thus the fold angle is not constant along the folds, particularly for the square pleats. Furthermore, the way in which the fold angles change along the length of a fold depends on which configuration the hypar folds into. Thus, manipulating the gap width (equivalently, the target fold angle) along the length of the pleats is a potential method for biasing folding.

To analyze changing behavior, each pleat was discretized into 30 segments and the fold angle of each segment was calculated from the node coordinates of the pleat’s two neighboring edges. The fold angle of fold $v11\u2212v21$ is illustrated in Fig. 9 as an example. The fold angle was computed using perpendicular vectors drawn from fold $v11\u2212v21$ to fold $v12\u2212v22$ and fold $v10\u2212v20$ in the flat state. The displacement of these vectors were tracked as the hypar self-folds in simulation, and the final fold angle was calculated from the final displacement of the vectors. On either end of the fold, there are two ways of computing the fold angle because the shorter edge (fold $v12\u2212v22$) tapers into two trapezoidal legs. These legs curve inward or outward depending on the hypar’s final configuration. The fold angle is calculated by terminating the perpendicular vectors at the legs; this methodology adds a discontinuity at vertices $v12$ and $v22$. To achieve a continuous plot, fold $v12\u2212v22$ is extended and the fold angle from the fold extension rather than the trapezoid legs is calculated. This changing fold angle behavior may be harnessed to control the hypar’s configuration, as the increasing or decreasing behavior changes from one configuration to another.

To implement this biasing technique, we assign a different fold angle to each discrete segment on the pleat folds, with increasing gap sizes in the direction of the increasing fold angle for the particular target configuration. Gap widths were designed to increase linearly from a 0.5 mm gap to a 1.5 mm gap depending on the prescribed configuration.

The results for pleat biasing show good performance both in simulation and physical experiments. For the 1–3 configuration, the hypars folded correctly for 82% of the simulated trials; for the 2–4 configuration, the hypars folded correctly for 92% of the trials. With a *p*-value of <0.0001 for both configurations, pleat biasing is statistically significantly different from the unbiased hypar, indicating that it is successful in biasing the hypar into one configuration or another. In physical tests, although the results were not statistically significant, the pleat biasing approach did seem to achieve a higher percentage of correct configurations ($14/20=70%$ success).

### 4.3 Changing Fold Angles of Different Folds (Diagonals Biasing).

For our third biasing method, we take inspiration from Ref. [19], who analyzed the configuration space of a pattern. For this study, the hypar has many folds, and it is not practical to map the full configuration space. We instead consider pairs of diagonal folds. For example, plotting the fold angles of fold $v20\u2212v21$ and fold $v10\u2212v11$ against each other gives us a projection of the configuration space, as shown in Fig. 10. The configuration space shows two distinct, well-separated configuration curves. The figure also shows the trajectories taken by hypar to which randomness has been applied in simulation. Although these hypar deviate from each other, most are tightly clustered around the theoretical trajectories. The trajectories that follow paths diagonally down the middle of the plot indicate hypars that initially folded into indistinct configurations, before snapping into a particular configuration.

Based on these observations, we then programed target fold angles that lie along the tangents to those curves. This analysis indicates that applying a force in the flat state along fold $v10\u2212v11$ that is 2.5 times stronger than the force along fold $v20\u2212v21$ will produce a force that is tangent to the 1–3 configuration trajectory. Similarly a force that is 2.5 times stronger along the opposite diagonals should push the pattern towards the 2–4 configuration.

We computed tangent vectors for each of the pleats in the hypar and the ratios of target fold angles. These vectors were applied to the standard hypar crease pattern by modifying the equilibrium fold angle of the diagonal creases according to the model in Sec. 3.1. Thus, by scaling all equilibrium fold angles along the *v*_{1} − *v*_{3} diagonal by the corresponding scale factor, the hypar is biased to fold into the 1–3 configuration. Similarly, by scaling all equilibrium fold angles along the *v*_{2} − *v*_{4} diagonal, the hypar is biased to fold into the 2–4 configuration.

The results indicate that diagonal biasing is the most successful biasing method in simulation and achieves similar results to pleat biasing in physical experiments ($14/20=70%$ success). The simulation trials showed 100% success for both configuration, with no hypars failing to fold into a distinct configuration.

### 4.4 Discussion

#### 4.4.1 External Forces Present in Physical Tests.

Our results show that it is possible to bias a self-folding pattern to change the distribution of resulting folded states, and we have demonstrated three different approaches to doing this. While the pleat and diagonals biasing approaches did seem to have some positive effect on the reliability of self-folding, the center biasing approach failed in physical experiments. We also observe that in all of the methods, the physical tests showed a lower success rate than the simulation, despite our accounting for fabrication error in the simulation. This result is not entirely unexpected. Due to the compliance of the laminate, if any additional external forces act on the pattern, it is possible that the hypar will fold into the incorrect configuration anyway. In physical tests, the samples were placed in a water bath, which provided additional resistive forces on the folding structure that were not modeled in the simulation. Practically, during experiments, we discovered that the orientation that the pattern was dropped into the hot water bath was an important factor in the results. For all of these tests, the hypar was therefore lowered into the bath vertically to reduce as much as possible the effects of horizontal forces on the final configuration. The results show that in the future, it may also be useful to add random forces to the simulation to capture reliability of folding more accurately.

#### 4.4.2 Modeling the Effects of Contact and Delamination.

The layers in our self-folding structures are held together by an adhesive layer between the rigid PET layer and the heat-shrinking PVC. Practically, there will be delamination and shifting of layers as the structure self-folds, particularly when the rigid PET layers come into contact as the PVC layer shrinks. In simulations, we scaled the shrink ratio *λ* to compensate for the average effects of this problem and injected fold angle variation to address adhesion issues in individual patterns. The values used for scaling *λ* and for fold angle variation needed to be measured experimentally, however. In the future, it would be useful to characterize the effect of the adhesive layer and material collisions on the model, as well as to improve adhesion between layers in the fabrication process.

#### 4.4.3 Extensions to the Hypar Tessellation.

Finally, based on these results, we also tested whether the same approaches could be used to bias self-folding of a more complex pattern, a tessellation of hypars. The hypar tessellation is a composition of the single hypar patterns, so the purpose of this test was to see if biasing approaches could be also be “composed,” that is whether biasing methods that are successful on a small pattern could still be applied to the same pattern when it is a part of a larger structure. We chose the diagonals biasing method, which seemed to be the most successful in simulations and experiments. Using these patterns, we programed two different 2 × 2 tessellations, shown in Fig. 11. In one configuration, the units were arranged with alternating assignments, while in another, they were all arranged with the same assignment. Each design was fabricated and tested three times. The results for both designs was that two of the three samples folded correctly into the desired configurations, while one did not. For the alternating assignment, 11 of the 12 total units folded into the correct configuration, while for the constant assignment, 10 of the 12 units folded into the correct configuration. Similarly to our single hypar results, these results indicate great promise that fold biasing methods can be used in the future to improve folding reliability.

One additional observation is that for the tessellation tested here, although the configurations of the individual units were prescribed, each of the tessellations actually had two potential forms [26] since the orientation of the center vertex can also be controlled independently of the units. Incorporating vertex biasing methods such as those explored in Ref. [23] may be useful for this pattern in the future.

## 5 Conclusions

In this paper, we present a simulation approach to self-folding origami patterns and investigate the effect of various fold pattern modifications on the folded state of an origami hyperbolic paraboloid. We characterize the fold angle error in our self-folding approach and propose a model for fold angle and stiffness based on the gap width. Three methods of biasing the fold pattern are then tested, showing that the biasing methods could be used to improve folding reliability in simulation. The results include the effect of random noise from fabrication error. We also verify these findings in physical experiments, and the results show that the biasing methods may have an impact (although the results were not statistically significant), potentially increasing folding accuracy from 50% (purely random) to 70% (Table 4). Finally, we extend the results to a 2 × 2 hypar tessellation, showing the biasing methods could be used in more complex patterns as well. Future work includes incorporating the effects of environmental forces into our simulation, testing our biasing approaches on other patterns, and comparing these methods to other biasing approaches that include vertex manipulation to determine how to fold these patterns more robustly.

## Acknowledgment

Support for this project has been provided in part by NSF Grant No. 1659190. We also thank the Research Science Institute (RSI) program, including Dr. Amy Sillman, River Grace, and the RSI staff, for helpful feedback, discussions, and administrative support.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

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